 	program testefcn

        include 'param_dim.inc'

      	integer m,numpar,ldfjac	
	integer maxfev,mode,nprint,info,nfev,njev
	integer ipvt(npar_max)
	double precision ftol,xtol,gtol,factor

	double precision p(npar_max),s(npar_max)
	double precision x(npar_max),fvec(NRCX),fjac(NRCX,npar_max),
     *        diag(npar_max), qtf(npar_max), wa1(npar_max),
     *        wa2(npar_max), wa3(npar_max), wa4(NRCX)
	double precision epsfcn

        double precision fvec1(NRCX)

	double precision fjacDD(NRCX*npar_max),fjacDF(NRCX*npar_max)
	double precision errodif(npar_max),normadif


        double precision r,dudp(Nparg_max),h
	integer count1,count2,count_rate,count_max

	integer FDAT,FDCC,FDEX,FRES,FROT
	CHARACTER*15 FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES


	common/geom/p,s,numpc
	COMMON /narq/ FDAT,FDCC,FDEX,FRES,FROT
	COMMON /cFDAT/ FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES
! 	external fcn
	external fcn1,fcn2


	call input_opt(ftol,xtol,gtol,maxfev,numpar,
     &           m,metodo,mode,nprint,factor,epsfcn)
	ldfjac = NRCX


	if(numpar .eq. numpc)then
            do i=1,numpar
                !x(i)=s(i)
                x(i)=tan((s(i)-0.5)*3.1416)
            enddo
	else if(numpar .eq. 2*numpc)then
            do i=1,numpar
                x(i)=p(i)
            enddo
	else if(numpar .eq. 3*numpc)then
            do i=1,numpc
                x(3*i-2)=p(2*i-1)
                x(3*i-1)=p( 2*i )
                !x(3*i  )=s(  i  )
                x(3*i  )=tan((s(i)-0.5)*3.1416)
            enddo
	endif





! ! Analise da influencia do epsfcn na norma do vetor (fvec1-fvec)/h
! 
!         iflag = 1
!         metodo = 1
!         
!         r = RESID(x,fvec,fjac,m,numpar,m,iflag,metodo)
!         
! 	iflag = 2
!         metodo = 2
!         
!         r = RESID(x,fvec,fjacDD,m,numpar,m,iflag,metodo)       
! 
! 	open(1,file='hj-norma-jacDF.txt')
! 	open(2,file='hj-norma-erro-jac.txt')
!         iflag = 1
!         metodo = 1        
!         do k=0,28
! 
!             epsfcn = 2.**(-12.+k*.5)
! 
!             h = dsqrt(epsfcn)
!             
!             do j = 1, numpar
!                 temp = x(j)
!                 
!                 x(j) = temp + h
!                 r = RESID(x,fvec1,fjac,m,numpar,m,iflag,metodo)
!                 x(j) = temp
!                 
!                 dudp(j) = 0.
! 		errodif(j) = 0.
! 		normadif = 0.
!                 do i = 1, m
!                     dudp(j) = dudp(j) + ((fvec1(i) - fvec(i))/h)**2
! 		    fjacDF((j-1)*m+i) = (fvec1(i) - fvec(i))/h	
! 		    errodif(j) = errodif(j) + (fjacDD((j-1)*m+i) - 
!      &					fjacDF((j-1)*m+i))**2
! 		    normadif = normadif + fjacDD((j-1)*m+i)**2
!                 enddo
!                 dudp(j) = dsqrt(dudp(j))
!     		errodif(j) = dsqrt(errodif(j))/normadif
!             enddo
!         
!             write(1,6)epsfcn,(dudp(j),j=1,numpar)
! 	    write(2,6)epsfcn,(errodif(j),j=1,numpar)
!         enddo
! 	close(1)
! 	close(2)
! 
!  6      format(f14.8,12f12.6)
!         stop
! !   fim da analise da influencia do epsfcn no calculo do jacobiano



	if(metodo.eq.1)then

            write(4,*)'_______________________________________'
	    write(4,*)'--L-M - J:Finite Differences--'
            write(4,*)'_______________________________________'

            call system_clock(count1,count_rate,count_max)

	    call lmdifmodif(fcn1,m,numpar,x,fvec,ftol,xtol,gtol,maxfev,
     *            epsfcn,diag,mode,factor,nprint,info,nfev,fjac,ldfjac,
     *                ipvt,qtf,wa1,wa2,wa3,wa4)
            call system_clock(count2,count_rate,count_max)
            r=(count2-count1)/count_rate

            write(4,*)'Fim da primeira etapa da otimizacao'
            write(4,*)'Tempo de resol. (1a etapa) c/ LMDIF:',r,' s'
            write(4,*)'info: ',info,'  nfev: ',nfev
	else if(metodo.eq.2)then

            write(4,*)'_______________________________________'
            write(4,*)'--L-M - J:DirectDifferentiation--'
            write(4,*)'_______________________________________'

            call system_clock(count1,count_rate,count_max)

            call lmder(fcn2,m,numpar,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
     *                 maxfev,diag,mode,factor,nprint,info,nfev,njev,
     *                 ipvt,qtf,wa1,wa2,wa3,wa4)
            call system_clock(count2,count_rate,count_max)
            r=(count2-count1)/count_rate

            write(4,*)'Tempo de resol. (1a etapa) c/ LMDER:',r,' s'
            write(4,*)'info: ',info,'  nfev: ',nfev,'  njev:',njev
	endif




        !  2o passo da otimizacao   - 
	if(numpar .eq. 2*numpc)then
            write(4,*)'************************************'
            write(4,*)' 2o passo da otimizacao - x y s '
            write(4,*)'************************************'

!             xtol = xtol/2.
!             ftol = ftol/2.
!             gtol = gtol/2.

            do i=1,numpar
                p(i)=x(i)
            enddo
            
	
            numpar = 3*numpc

		write(4,*)
		write(4,*)'  Parametros do alg. de otimizacao'
		write(4,*)
		write(4,*)'   ftol: ',ftol
		write(4,*)'   xtol: ',xtol
		write(4,*)'   gtol: ',gtol
		write(4,*)'   maxfev: ',maxfev
		write(4,*)'   numpar: ',numpar		

            do i=1,numpc
                x(3*i-2)=p(2*i-1)
                x(3*i-1)=p( 2*i )
                x(3*i  )=tan((s(i)-0.5)*3.1416)
            enddo
    
            if(metodo.eq.1)then

                write(4,*)'_______________________________________'
                write(4,*)'--L-M - J:Finite Differences--'
                write(4,*)'_______________________________________'


                call system_clock(count1,count_rate,count_max)

                call lmdifmodif(fcn1,m,numpar,x,fvec,ftol,xtol,gtol,
     *                   maxfev,epsfcn,diag,mode,factor,nprint,info,
     *                   nfev,fjac,ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
                call system_clock(count2,count_rate,count_max)
                r=(count2-count1)/count_rate

                write(4,*)'Tempo de resol. (2a etapa) c/ LMDIF:',r,' s'
                write(4,*)'info: ',info,'  nfev: ',nfev

            else if(metodo.eq.2)then

                write(4,*)'_______________________________________'
                write(4,*)'--L-M - J:DirectDifferentiation--'
                write(4,*)'_______________________________________'   

                call system_clock(count1,count_rate,count_max)

                call lmder(fcn2,m,numpar,x,fvec,fjac,ldfjac,ftol,xtol,
     *                 gtol,maxfev,diag,mode,factor,nprint,info,nfev,
     *                 njev,ipvt,qtf,wa1,wa2,wa3,wa4)
                call system_clock(count2,count_rate,count_max)
                r=(count2-count1)/count_rate

                write(4,*)'Tempo de resol. (2a etapa)  c/ LMDER:',r,' s'
                write(4,*)'info: ',info,'  nfev: ',nfev,'  njev:',njev

            endif	

        endif !fim da segunda parte da otimizacao

 5	format(i5,54f10.5)

	close(4) ! arquivo xxxx.dop.out
	close(FRES)
	
	
	open(1,file='snr.tmp')
	read(1,*)snr
	close(1)
	
	open(1,file='pgeom.tmp')
	read(1,*)g1,g2,g3,g4,g5,g6,g7
	close(1)
	
	open(1,file='fobj0.tmp')
	read(1,*)fobj0
	close(1)
	
	open(1,file='fobjf.tmp')
	read(1,*)fobjf
	close(1)
	
	open(1,file='resum-'//FILE_OPT//'.txt')
	write(1,*)snr,g1,g2,g3,g4,g5,g6,g7,fobj0,fobjf
	close(1)	
	
	call system('gzip -f '//FILE_RES)
	call system('rm *.tmp ')

	stop
	end

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! ! !  PARA USAR COM LMDIF - versao secante do metodo de levenberg marquardt
	subroutine fcn1(m,n,x,fvec,iflag)
	include 'param_dim.inc'
	double precision x(n)
	double precision fobj
	integer iavalfunct,n,m,ldfjac,metodo
	double precision fvec(m),fjac(NRCX*npar_max)
        integer count1,count2,count_rate,count_max
	data iavalfunct/0/
	iavalfunct=iavalfunct+1


	metodo =1
	if(iflag .eq. 0)then

            write(4,*) '_______________________________________'
            fobj=resid(x,fvec,fjac,m,n,ldfjac,iflag,metodo)    
            write(4,*) '_______________________________________'

	elseif(iflag .eq. 1)then

            call system_clock(count1,count_rate,count_max)
!             write(4,*) 'fcn1 - calculo do vetor residuo - FVEC'
            fobj=resid(x,fvec,fjac,m,n,ldfjac,iflag,metodo)    
            call system_clock(count2,count_rate,count_max)
            r=(count2-count1)
            write(4,*)'tempo avaliacao: (s)',r/count_rate
            write(4,*)
	elseif(iflag .eq. 2)then

            call system_clock(count1,count_rate,count_max)
!             write(4,*) 'fcn1 - calculo do jacobiano de f - FJAC'
            fobj=resid(x,fvec,fjac,m,n,ldfjac,iflag,metodo)    
            call system_clock(count2,count_rate,count_max)
            r=(count2-count1)
            write(4,*)'tempo aval p grad: (s)',r/count_rate
            write(4,*)
	endif



! 	write(*,*)'x: '
! 	write(*,1000)(x(i),i=1,n)
! 	write(*,*)'fvec'
! 	write(*,1001)(fvec(i),i=1,m)

 1000	format(8f10.4)
 1001	format(104f10.4)

 1130   format(22x,3f14.6)
	return
	end	
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! ! !  PARA USAR COM LMDER
	subroutine fcn2(m,n,x,fvec,fjac,ldfjac,iflag)
	include 'param_dim.inc'
	real*8 x(n),xplus(npar_max)
	double precision fobj
	integer iavalfunct,n,m,ldfjac,metodo
	double precision fvec(m),fjac(ldfjac,n),fvecplus(NRCX)

	real*8 fvecaux(NX),xaux(npar_max)

        integer count1,count2,count_rate,count_max

	data iavalfunct/0/
	iavalfunct=iavalfunct+1

	metodo = 2
	if(iflag .eq. 0)then

            write(4,*) '_______________________________________'
            fobj=resid(x,fvec,fjac,m,n,ldfjac,iflag,metodo)
            write(4,*) '_______________________________________'

	else if(iflag .eq. 1)then
            call system_clock(count1,count_rate,count_max)
!             write(4,*) 'fcn2 - calculo do vetor residuo - FVEC'
            fobj=resid(x,fvec,fjac,m,n,ldfjac,iflag,metodo)
            call system_clock(count2,count_rate,count_max)
            r=(count2-count1)
            write(4,*)'tempo avaliacao: (s)',r/count_rate
            write(4,*)

	else if(iflag .eq. 2)then
            call system_clock(count1,count_rate,count_max)
!             write(4,*) 'fcn2 - calculo do jacobiano de f - FJAC'
            fobj=resid(x,fvec,fjac,m,n,ldfjac,iflag,metodo)
            call system_clock(count2,count_rate,count_max)
            r=(count2-count1)
            write(4,*)'tempo Jacobiano: (s)',r/count_rate
            write(4,*)
	endif

 1000	format(8f10.4)
 1001	format(104f10.4)
 300    format(10f10.6)
 1130   format(22x,3f14.6)
	return
	end	

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 

C
C funcao para calculo de problemas de potencial bi-dimensionais pelo
C metodo dos elementos de contorno com elementos constantes.
C

	FUNCTION RESID(PAR,fvec,dudQOT,m,numpar,ldfjac,iflag,metodo)
	
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)


	INCLUDE 'param_dim.inc'

	
	CHARACTER*11 FILE_B
	CHARACTER*15 FILE_DAT,FILE_DCC,FILE_DEX,FILE_RES,
     &             FILE_ROT,FILE_OPT
	character*19 FILE_OPTR
	INTEGER FDAT,FDCC,FDEX,FRES,FROT,FOPTR
	integer iterp,direcao,rotina
	real resid
	double precision p(npar_max),s(npar_max)

	DIMENSION GF(3*NX/5,3*NX/5),HF(3*NX/5,3*NX/5)
	DIMENSION XF(3*NX/5),YF(3*NX/5),INCF(3*NX/5,2)
	DIMENSION G(NX,NX),H(NX,NX),Hzero(NX,NX)
	DIMENSION X(NX),Y(NX),XM(NX),YM(NX),FI(NX),DFI(NX)
	DIMENSION INC(NX,2),KODE(NX),IT(NX)
	DIMENSION RCONT(NRCX),RCEX(NRCX),fvec(m)

	DIMENSION XIS(NX),U(NX,NBX+1),HBAR(NBX+1,NBX+1),EBAR(NBX+1)
	DIMENSION EEBAR(NBX+1),YIP(NBX+1),CE(NBX+1),ES(NBX+1),XISA(NX)

	DIMENSION DB(NX,NC2MX),NP(NX)
	INTEGER IDCASO(NCasosX,3),I_ELET(16,NELEX)

	DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max), NUMPCI(ninc_max)
	integer VARIA(ninc_max)

	double precision dudQ(NX*Nparg_max)
        dimension DFIaux(NX)
	double precision dudQOT(ldfjac,numpar)

	real*8 PAR(numpar),PARG(npar_max)

	integer PCnnp(nparg_max) ! conta qtos nohs parametros cada ponto de controle influencia
	integer PCnp(nparg_max,NX) ! guarda quais nohs parametros cada ponto de controle influencia
	double precision PCgama(nparg_max,NX) ! guarda os pesos com que cada coord dos PC influencia cada noh parametro

	double precision PCgamaSx(nparg_max,NX),PCgamaSy(nparg_max,NX)! guarda os pesos com que cada parametro S dos PCs influencia cada Coordenada X ou Y dos nohs parametros

	integer jNP,cjNPx,cjNPy
	double precision DB2(NX,NC2MX)
	double precision DFI0(NX),DFI1(NX),FI1(NX)
	DIMENSION G1(NX,NX),H1(NX,NX)
	logical viavel
	double precision fvecsempena(NRCX)

        real tarray(2),t1sist,t2sist,t1grad,t2grad


        ! Para a rotina de SVD
	real*8 A(NX,NX),Sigma(NX),WORK(5*NX),epssvd
	real*8 Usvd(NX,NX),VT(NX,NX),Vsvd(NX,NX)
	CHARACTER*1 JOBU,JOBV
        real*8 A2(NX,NX),Ainv(NX,NX)

	! para a verificacao da sol do sistema
	real*8 bsist(NX),vRsist(NX)


        double precision Haux(NX,NX),Eaux(NX)
        integer IPIV(NX)

	COMMON /narq/ FDAT,FDCC,FDEX,FRES,FROT
	COMMON /cFDAT/ FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES
   
	COMMON /PG/ GI(8,2),OME(8,2)
	common/plota/ x,y,inc,NN,nfixaux,i_elet,nepe
	common/idr/iterp,direcao,rotina
	
	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA
	common/geom/p,s,numpc
	common/dadossensib/PCnnp,PCnp,PCgama,PCgamaSx,PCgamaSy


	save iaval,GF,HF,XF,YF,XM,YM,RCEX,FI,DFI,INCF,NFIX
	save VCC,NCasos,IDCASO,ITCC,NG
	save fvecsempena,inviacont

	integer count1,count2,count_rate,count_max,r



	DATA
     *   (GI(J,1),J=1,4)/0.86113631,-0.86113631,0.33998104,-0.33998104/,
     *   (OME(J,1),J=1,4)/0.34785485,0.34785485,0.65214515,0.65214515/,
     *   (GI(J,2),J=1,8)/0.96028986,-0.96028986,0.79666648,-0.79666648,
     *                   0.52553210,-0.52553210,0.18343464,-0.18343464/,
     *   (OME(J,2),J=1,8)/0.10122854,0.10122854,0.22238103,0.22238103,
     *                    0.31370665,0.31370665,0.36268378,0.36268378/

	data iaval/1/
	data iplota/0/

	! Constante alterando a solucao fundamental
	ALFA = 100.0 
	
	PI= 4.D0*ATAN(1.D0)
	B=  0.D0
	xl_el = .5

	! Parametros para solucao do sistema.

	ITRI= 0
	! WRITE(FRES,110) 
 110  	FORMAT(//,80('*'),/,80('*'),//)


	! Informacoes sobre o numero da entrada de dados,e o numero da
	! saida de resultados.

	FDAT= 12
	FDCC= 13
	FDEX= 14
	FRES= 16
	FROT= 17
	FOPTR=19


	! Monta nomes do arquivo de saida.    
 
	FILE_B   = FILE_DAT
	FILE_DCC = FILE_B//'.dcc'
	
	FILE_B   = FILE_OPT 
	FILE_RES = FILE_B//'.res'
	
	FILE_B   =FILE_OPT
	FILE_OPTR=FILE_B//'_OPT.res'
	
	
	OPEN (FDAT,FILE=FILE_DAT,STATUS='OLD')
	OPEN (FDCC,FILE=FILE_DCC,STATUS='OLD')
	OPEN (FDEX,FILE=FILE_DEX,STATUS='OLD')
	OPEN (FRES,FILE=FILE_RES,STATUS='UNKNOWN')
	open (FOPTR,FILE=FILE_OPTR,status='unknown')


	! INPUT
	IF(IAVAL.EQ.1) THEN
            ! Le os dados globais para as analises
            FILE_ROT = FILE_B//'.rot'
            OPEN (FROT,FILE=FILE_ROT,STATUS='UNKNOWN')

            CALL INPUT0(X,Y,RCEX,VCC,INC,NFIX,NEPE,I_ELET,
     &              NCasos,IDCASO,ITCC,NG,NE)

 125  	    format(5x,'Parametro da solucao fundamental:',e12.5)

            DO I=1,NFIX
                XF(I)= X(I)
                YF(I)= Y(I)
                INCF(I,1)= INC(I,1)
                INCF(I,2)= INC(I,2)
            ENDDO
	ELSE
            DO I=1,NFIX
                X(I)= XF(I)
                Y(I)= YF(I)
                INC(I,1)= INCF(I,1)
                INC(I,2)= INCF(I,2)
            ENDDO
	ENDIF


	if(numpar .eq. numpc)then
            do i=1,numpc
                parg(i*3-2) = p(i*2-1)
                parg(i*3-1) = p(i*2  )
                parg(i*3  ) = (atan(par(i)))/3.1416+0.5
            enddo
	else if(numpar .eq. 2*numpc)then
            do i=1,numpc
                parg(i*3-2) = par(i*2-1)
                parg(i*3-1) = par(i*2  )
                parg(i*3  ) = s(i)
            enddo
	else if(numpar .eq. 3*numpc)then
            do i=1,numpc
                parg(i*3-2) = par(i*3-2)
                parg(i*3-1) = par(i*3-1)
                parg(i*3  ) = (atan(par(i*3)))/3.1416+0.5
            enddo
	endif



	if(iflag .eq. 1 )then
		viavel = .true.

! 		call viabilidade(PARG,X,Y,numpar,numpc,nfix,viavel,xl_el)

                if(iaval .eq. 1 .and. .not. viavel)then
			write(*,*)'Aproximacao inicial nao viavel'
			write(4,*)'Aproximacao inicial nao viavel'
			stop
		endif


		if(viavel)then
			pena = 1.
		else

			write(4,*)'INVIAVEL pena: ',1.+pena
			write(*,*)'INVIAVEL pena: ',1.+pena
	
			do i=1,m
			fvec(i)=fvecsempena(i)*(1.+pena) !penalizacao
			enddo
			iaval = iaval+1
			pena = pena/2.
			return
	    	endif
	endif

!         write(*,*)'parg: ',(PARG(i),i=1,numpc*3)
	CALL GERACONT(PARG,XL_EL,X,Y,NUMPC,FRES,INC,NFIX,NN,N_eq,
     &					npa0) 

	nfixaux = NFIX
	write(frot,*) '      Avaliacao:',iaval


        if(iflag.eq.0)then !só plotagem
 	    write(4,*) '                PLOTAGEM No. ',iplota
            resid = 0.
            do i=1,m
                resid = resid + fvec(i)**2
            enddo
            resid = resid/2.0
                
            write(4,222)(PARG(i),i=1,3*numpc)
            write(4,*) 'Residuo:',RESID
            write(4,*)


            call plota_cont(PARG,resid,iplota,numpar,numpc)
            iplota = iplota +1
            return
        endif

        DO IC= 1,NCasos
            ! Introduz as condicoes de contorno com base no caso de carregamento
            WRITE(FRES,1110) 
 1110       FORMAT(//,80('*'),/,80('*'),//)
            WRITE(FRES,1120) IC
 1120       FORMAT(5X,'CASO DE CARREGAMENTO: ',I3)

            CALL DETCC(FI,DFI,VCC,FRES,NEPE,I_ELET,IC,IDCASO,NFIX,KODE)

            ! Formar o sistema de equacoes.
            CALL FMAT(G,H,GF,HF,X,Y,XM,YM,FI,DFI,B,alfa,
     &                KODE,INC,NG,nfix,NN,N_eq,iaval)


!            !!Solucoes do sistema de equacoes.
! 
!                 ! 1) Eliminacao de gauss
!                 call system_clock(count1,count_rate,count_max)
!                 ITRI=0
!                 nc2m = 1
!                 CALL SLNPDM(H,DFI,etol,det,IT,itri,N_eq,nc2m)
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'SLNPDM(s):',r,'/rate',' NRHS: ',
!      &                  1,' Eqs: ',N_eq 
! 		raux=0.
! 		do i=1,N_eq
! 			raux=raux+dfi(i)**2
! 		enddo
! 		write(*,*)'norma sol: ',dsqrt(raux)
! 		stop
!                 ! fim do solve com eliminacao de gauss


!                 ! 2) GMRES
!                 call system_clock(count1,count_rate,count_max)
!                 NB=    NBX
!                 LMAX=  3
!                 ICONV= 0
!                 IRP=   0
!                 CALL GMRES(H,DFI,XIS,XISA,U,EBAR,EEBAR,HBAR,YIP,CE,
!      &            ES,ETOL,N_eq,NB,NB,LMAX,ICONV,ITER,IRP,FRES)
!                 IF(ICONV.EQ.0) STOP'CONVERGENCIA NAO ATINGIDA - gmres'
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'GMRES(s):',r,'/rate',' NRHS: ',
!      &                  1,' Eqs: ',N_eq 
!                 !fim do solve com GMRES



                !3) Solve com subrotina do lapack
!                 call system_clock(count1,count_rate,count_max)

                NRHS = 1
                call DGESV(N_eq,NRHS,H,N_eq,IPIV,DFI,N_eq,INFOsolv)

! 			raux=0.
! 			do i=1,N_eq
! 				raux=raux+dfi(i)**2
! 			enddo
! 			write(*,*)'norma sol: ',dsqrt(raux)
!  
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'DGESV(s):',r,'/rate',' NRHS: ',
!      &                1,' Eqs: ',N_eq 

               ! fim do solve com lapack



C! ! ! ! ! ! 		! processo iterativo para melhor a sol de sistemas com matrizes mal cond (lu- lapack)(GOLUB)
C! ! ! ! ! !		! lembrar de fazer copia da matriz e do vetor dfi antes de resolver o sistema 
C! ! ! ! ! ! 		nvezes = 10
C! ! ! ! ! ! 		do k=1,nvezes
C! ! ! ! ! ! 			
C! ! ! ! ! ! 			call matvec(A,DFI,vRsist,N_eq)
C! ! ! ! ! ! 			do i=1,N_eq
C! ! ! ! ! ! 				vRsist(i) = bsist(i) - vRsist(i)
C! ! ! ! ! ! 			enddo
C! ! ! ! ! ! 
C! ! ! ! ! ! 			raux=0.
C! ! ! ! ! ! 			do i=1,N_eq
C! ! ! ! ! ! 				raux=raux+vRsist(i)**2
C! ! ! ! ! ! 			enddo
C! ! ! ! ! ! 			write(*,*)'    norma do resid da sol: ',dsqrt(raux)
C! ! ! ! ! ! 
C! ! ! ! ! ! 			CALL DGETRS( 'No transpose', N_eq,NRHS,H,N_eq,IPIV, 
C! ! ! ! ! !      &                			vRsist,N_eq,INFOsolve )
C! ! ! ! ! ! 
C! ! ! ! ! ! 			do i=1,N_Eq
C! ! ! ! ! ! 				dfi(i) = dfi(i) + vRsist(i)
C! ! ! ! ! ! 			enddo 
C! ! ! ! ! ! 
C! ! ! ! ! ! 			raux=0.
C! ! ! ! ! ! 			do i=1,N_eq
C! ! ! ! ! ! 				raux=raux+dfi(i)**2
C! ! ! ! ! ! 			enddo
C! ! ! ! ! ! 			write(*,*)'norma sol: ',dsqrt(raux)
C! ! ! ! ! ! 
C! ! ! ! ! ! 		enddo
C! ! ! ! ! ! 		! fim do processo que melhoraria a solucao(Golub)



!              ! 4) Solve com a pseudoinversa para teste:
!                ! verifica a matriz de coeficientes com svd
! 
!                 do i=1,NX
!                     do j=1,NX
!                         A(i,j)=H(i,j)
!                     enddo
!                enddo
! 
! 		epssvd = 0.000001 !valores de sigma da SVD abaix desta tol sao tratados como zero
! 
! 		CALL MPINV (N_eq,N_eq,N_eq,N_eq,A,AINV,
!      &			Sigma,Eaux,Usvd,Vsvd,WORK,IERR,epssvd)
! 
!                 call solvepseudo(Ainv,dfi,N_eq,1)
! 
! 		write(*,*)'Sigma: ',(Sigma(i),i=1,N_eq)
! 		stop
!                ! fim do solve com pseudoinversa da SVD

 
            if (metodo.eq.2 .and. iflag.eq.2)then



!                 call system_clock(count1,count_rate,count_max)


                NNP = NN - NFIX !todos os nohs das splines internas sao nohs parametro
			
                do inp=1,NNP
                        NP(inp) = NFIX + inp ! numeros dos nohs parametros
                enddo

! 		Formar o 2o membro do problema das sensibilidades.
                CALL FMATD(X,Y,XM,YM,DB,FI,DFI,alfa,KODE,INC,NG,NNP,NP,
     &						NFIX,NN,N_eq)


                call calcDuDQ(dudQ,DB,PAR,numpc,numpar,N_eq,NFIX,NNP)


!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'montagem(s):',r,'/rate'





!                 !Solve com a pseudoinversa para teste:
!                 call solvepseudo(Ainv,dudq,N_eq,numpar)



!                 ! 1) Eliminacao de gauss
!                 call system_clock(count1,count_rate,count_max)
! 
!                 ITRI=1 !matriz jah triangularizada
!                 nc2m = numpar ! numero de sistemas a resolver
!                 CALL SLNPDM(H,dudQ,etol,det,IT,itri,N_eq,nc2m)
! 
!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'SLNPDM(s):',r,'/rate',' NRHS: ',
!      &                  numpar,' Eqs: ',N_eq 
! 
!                 ! fim do solve com eliminacao de gauss


!                 !2) GMRES
!                 call system_clock(count1,count_rate,count_max)

!                 call preGMRES(H,dudQ,N_eq,numpar,ITER,FRES)

!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'GMRES(s):',r,'/rate',' NRHS: ',
!      &                      numpar,' Eqs: ',N_eq 
!                 ! fim do solve com GMRES    



                !3) Solve com subrotina do lapack - para matria H já triangularizada
!                 call system_clock(count1,count_rate,count_max)

                NRHS = numpar
                CALL DGETRS( 'No transpose', N_eq,NRHS,H,N_eq,IPIV, 
     &                dudQ,N_eq,INFOsolve )

!                 call system_clock(count2,count_rate,count_max)
!                 r = (count2-count1)
!                 write(*,*)'DGESV(s):',r,'/rate',' NRHS: ',
!      &                  NRHS,' Eqs: ',N_eq 

                !fim do solve com lapack




                call sensib_OUTPT_OT(dudQ,dudQOT,numpar,numpc,
     &               ldfjac,I_ELET,IDCASO,IC,NEPE,NCasos,N_eq)


!                 write(*,*)'tempo p calc. do grad.:',(t2grad-t1grad)

            endif !iflag =2

            ! SAIDA DE DADOS
            CALL OUTPT(XM,YM,FI,DFI,KODE,NFIX,NN,N_eq)	

            ! Saida de dados para a otimizacao.
            CALL OUTPT_OT(FI,DFI,RCONT,RCEX,RESID,I_ELET,
     &                IDCASO,IC,NEPE,NCasos,FDEX,FROT)
        ENDDO

! 	WRITE(*,1010) ((RCONT(I)),I=1,NRC)
 1010   FORMAT(13E12.4,/,7(13E12.4,/))

	
!  calculo do vetor residuo
        if(metodo.eq.1  .or. (metodo.eq.2 .and. iflag.eq.1))then
            do i=1,m
                fvec(i)=RCONT(i)-RCEX(i)
                fvecsempena(i) = fvec(i)
            enddo
        endif


!  saida 
	if(iflag.eq.1 .or.iflag .eq.2)then
	    write(4,*) '                        Avaliacao:',iaval
            if(numpar .eq. 2*numpc)then
                write(4,222)(PARG(i),i=1,3*numpc)
            else if(numpar .eq. 3*numpc )then
                do i=1,numpc
                    write(4,224)parg(i*3-2),parg(i*3-1),parg(i*3),
     &                                  par(i*3)
                enddo
            elseif(numpar .eq. numpc)then
                do i=1,numpc
                    write(4,224)parg(i*3-2),parg(i*3-1),parg(i*3),
     &                                  par(i)
                enddo
            endif
            write(4,*) 'Residuo:',RESID

            if(iaval.eq.1) write(foptr,1999)
            iaval = iaval+1
            write(FOPTR,*) iaval,resid,(PARG(i),i=1,3*numpc)
        endif
! 	
! 	call plota_cont(PARG,resid,iplota,numpar,numpc)
!  	stop

 5    format(54f10.5)
 222  format('                      ',3f10.4)
 224  format('                      ',3f10.4,f20.4)
 1999 format('#   Aval.    Resid.       It.Pow.   Direcao   Rotina   ',
     &       ' Num.Eq.   It.GMRES')
 2000 format(i5,5x,f10.6,5i10)

      CLOSE (FDAT)
      CLOSE (FDCC)
      CLOSE (FDEX)
!       CLOSE (FRES)


      RETURN
      END



        subroutine solvepseudo(Ainv,dudq,N_eq,numpar)

        real*8 Ainv(N_eq,N_eq)
        real*8 dudq(N_eq,numpar)
        real*8 sol(1000,100)


        do i=1,N_eq
            do j=1,numpar
                sol(i,j) = 0.
            enddo
        enddo        

        do k=1,numpar
            do i=1,N_eq
                do j=1,N_eq
                    sol(i,k) = sol(i,k) +  Ainv(i,j)*dudq(j,k)
                enddo
            enddo
        enddo

!         do k=1,numpar
!             r = 0.
!             do i=1,N_eq
!                 r = r + sol(i,k)**2
!             enddo
! !             write(*,*)'(x = A-1 * b) norma de x=  ',r
!         enddo

        do k=1,numpar
            do i=1,N_eq
                dudq(i,k)=sol(i,k)
            enddo
        enddo	

        return
        end

        subroutine pseudoinverse(Sigma,U,VT,Aaux,N_eq)

        real*8 Sigma(N_eq)
        real*8 U(N_eq,N_eq),VT(N_eq,N_eq),VTT(1000,1000)
        real*8 Aaux(N_eq,N_eq)
!         real*8 APaux(1000,1000)


        ! Verificar se A = U * Sigma * V^T:

!         r = 0.
!         do i=1,N_eq
!             do j=1,N_eq
!                 Aaux(i,j) = 0.
!                 APaux(i,j) = 0.
!                 r = r + A(i,j)**2
!             enddo
!         enddo

! !         write(*,*) 'Norma de A: ',r
! 
!         do i=1,N_eq
!             do j=1,N_eq
!                 if(i .ne. j)then
!                     Sigaux(i,j) = 0.0
!                 else
!                     Sigaux(i,j) = Sigma(i)
!                 endif
!             enddo
!         enddo
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     APaux(i,k) = APaux(i,k) + Sigaux(i,j)* VT(j,k) 
!                 enddo
!             enddo
!         enddo
! 
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     Aaux(i,k) = Aaux(i,k) + U(i,j)* APaux(j,k) 
!                 enddo
!             enddo
!         enddo

!         r = 0.
!         do i=1,N_eq
!             do j=1,N_eq
!                 r = r + Aaux(i,j)**2
!             enddo
!         enddo
! !         write(*,*) 'Norma de A2: ',r




        ! Calculo da pseudoinversa


!         do i=1,N_eq
!             do j=1,N_eq
!                 VTT(i,j) = VT(j,i)
!             enddo
!         enddo

	do j=1,N_eq
        	do i=1,N_eq
			if (abs(Sigma(j)) .lt. 0.0001)then
				VT(j,i) = 0.0
			else
				VT(j,i) = VT(j,i)*1.0/Sigma(j)
			endif
		enddo
        enddo




	do k=1,N_eq
		do i=1,N_eq
			Aaux(i,k) = 0.0
			do j=1,N_eq
				Aaux(i,k) = Aaux(i,k) + VT(j,i)* U(k,j) ! U(k,j) jah é a transposta de U(j,k)
			enddo
		enddo
	enddo







! 
!         do i=1,N_eq
!             do j=1,N_eq
!                 Aaux(i,j) = 0.
!                 APaux(i,j) = 0.
!             enddo
!         enddo
! 
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     if(i.ne.j)then
!                         APaux(i,k) = APaux(i,k) + Sigaux(i,j)* U(k,j)
!                     else
!                         if(Sigaux(i,i) .lt. 0.00001)then
!                             APaux(i,k) = APaux(i,k) + 0.0*U(k,j)
!                         else
!                             APaux(i,k) = APaux(i,k) + 
!      &                          1.0/Sigaux(i,j)*U(k,j)
!                         endif
!                     endif
!                 enddo
!             enddo
!         enddo
! 
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     Aaux(i,k) = Aaux(i,k) + VT(j,i)* APaux(j,k) !esta eh a pseudoinversa
!                 enddo
!             enddo
!         enddo


        return
        end




        subroutine preGMRES(H,dudQ,N_eq,numpar,ITER,FRES)       

        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        INCLUDE 'param_dim.inc'


        DIMENSION XIS(NX),U(NX,NBX+1),HBAR(NBX+1,NBX+1),EBAR(NBX+1)
        DIMENSION EEBAR(NBX+1),YIP(NBX+1),CE(NBX+1),ES(NBX+1),XISA(NX)	

        double precision H(N_eq,N_eq),dudq(N_eq,numpar),DFI(NX)

        integer fres,iter

! 
!         write(*,*)'Antes solve gmres'
!         do i=1,N_eq
!             write(*,*)(dudq(i,j),j=1,numpar)
!         enddo


!         do j=1,numpar
!             r=0.
!             do i=1,N_eq
!                     r = r + dudq(i,j)**2
!             enddo
!             write(*,*)j,'Antes solve gmres |b|= ',r
!         enddo



        do j=1,numpar

            do i=1,N_eq
                dfi(i) = dudQ(i,j)
            enddo

            !2) GMRES
            NB=    NBX
            LMAX=  3
            ICONV= 0
            IRP=   0
            CALL GMRES(H,dfi,XIS,XISA,U,EBAR,
     &            EEBAR,HBAR,YIP,CE,ES,ETOL,N_eq,NB,NB,LMAX,
     &              ICONV,ITER,IRP,FRES)
            IF(ICONV.EQ.0)STOP'CONVERGENCIA NAO ATINGIDA - gmres'
            !fim do solve com GMRES


            do i=1,N_eq
                dudQ(i,j) = dfi(i) 
            enddo

        enddo
    

! !         write(*,*)j,'Apos solve gmres'
!         do j=1,numpar
!             r=0.
!             do i=1,N_eq
!                     r = r + dudq(i,j)**2
!             enddo
! !             write(*,*)j,'Apos solve gmres |b|= ',r


!         enddo



        return
        end



	subroutine viabilidade(PARG,Xcont,Ycont,numpar,numpc,
     &				nfix,viavel,xl_el)

	include 'param_dim.inc'
	logical viavel
	logical intersect
	double precision Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,area
	double precision Xcont(NX),Ycont(NX)
	double precision xl_el,u,v

	double precision PARG(numpc*3),CINC(ninc_max),C0
        double precision x(nparg_max),y(nparg_max),xc,yc,areaT
	double precision Ju,Jv,Juv,compr
        double precision areaD,areaD2 
	integer NNINC(ninc_max)
	integer NUMPI(ninc_max), NUMPCI(ninc_max)
	integer VARIA(ninc_max)


	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA
	save iteste,areaD
	data iteste/0/


	iteste = iteste +1
	viavel = .true.


! Testa se a orientacao dos pontos de controle esta no sentido horario:
	if(iteste.eq. 1)then
		ipcacum = 0
		do inc=1,NINCL
			do i=1,NUMPCI(inc)
				x(i) = PARG((i+ipcacum)*3-2)
				y(i) = PARG((i+ipcacum)*3-1)
			enddo
		
! 			call centroide(x,y,NUMPCI(inc),0,xc,yc,areaT)

      			call propgeom(x,y,NUMPCI(inc),0,compr,areaT,
     &					xc,yc,Ju,Jv,Juv)
			areaT = - areaT


			ipcacum = ipcacum + NUMPCI(inc)
		
			if(areaT .gt. 0.0)then
				write(*,*)'Erro na orientacao dos pontos de controle'
				viavel = .false.
				return
			endif
		enddo
	endif





! Verifica se o contorno externo estah orientado corretamente
	if(iteste .eq. 1)then
! 		call centroide(Xcont,Ycont,nfix,0,xc,yc,areaD)

      		call propgeom(Xcont,Ycont,nfix,0,compr,areaD,
     &					xc,yc,Ju,Jv,Juv)
		areaD = - areaD

		if(areaD .lt. 0.)then !
			write(*,*)'Malha do contorno orientada incorretamente'
			stop
		endif
		areaD = abs(areaD)
	endif


! Verifica se todos os ptos de controle pertencem ao domínio
	do j=1,numpc
		areaD2 = 0.

		do i=1,nfix
			Ax = Xcont(i)
			Ay = Ycont(i)
			Bx = Xcont(i+1)
			By = Ycont(i+1)
			Cx = PARG(j*3-2)
			Cy = PARG(j*3-1)
			if(i.eq.nfix)then
				Bx = Xcont(1)
				By = Ycont(1)
			endif

			if(area(Ax,Ay,Bx,By,Cx,Cy).lt.0.)then
				write(*,*)'Ponto de controle fora do domínio!'
				viavel = .false.
				return
			endif
		enddo

	enddo




!! Testa se os segmentos de reta formados pelos pontos de controle nao se interceptam 

	ilatacum = 0
	do inc=1,NINCL
		
		ilat1 = ilatacum + 1
		ilat2 = ilatacum + NUMPCI(inc) 
		
		do ilat = ilat1 , ilat2
		
			ia = ilat
			if(ilat.eq.ilat2)then
				ib=ilat1
			else
				ib=ia+1
			endif
	
			Ax = PARG(ia*3-2)
			Ay = PARG(ia*3-1)
			Bx = PARG(ib*3-2)
			By = PARG(ib*3-1)
	! 	    

			jlatacum = 0
			do jinc = 1,NINCL 
				jlat1 = jlatacum + 1
				jlat2 = jlatacum + NUMPCI(jinc)
		    
				do jlat = jlat1,jlat2
					ic = jlat
					if(jlat.eq.jlat2)then
						id=jlat1
					else
						id=ic+1
					endif
			
					Cx = PARG(ic*3-2)
					Cy = PARG(ic*3-1)
					Dx = PARG(id*3-2)
					Dy = PARG(id*3-1)

					if(jlat .ge. ilat+2 .and. id .ne. ia)then

						if(intersect(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,u,v) .eqv. .true.)then
							viavel = .false.
							write(*,*)'1-INTERSECAO DOS SEGMENTOS ',ilat,jlat
							return
						elseif(u .le. xl_el .or. v .le. xl_el)then
							viavel = .false.
							write(*,*)'PC mto proximos',ilat,jlat
							return
						endif
					elseif(ib .eq. ic .and. jlat .gt. ilat)then

						if(theta(Ax,Ay,Bx,By,Dx,Dy,u,v) .lt. 3.14/18.)then ! 10graus entre as laterais
! 							viavel = .false.
! 							write(*,*)'Angulo theta<10graus ',ilat,jlat
! 							return
						elseif(u .le. xl_el .or. v .le. xl_el)then
							viavel = .false.
							write(*,*)'PC mto proximos',ilat,jlat
							return
						endif
					elseif(ia .eq. id .and. jlat .gt. ilat)then

						if(theta(Cx,Cy,Dx,Dy,Bx,By,u,v) .lt. 3.14/18.)then ! 10graus entre as
! 							viavel = .false.
! 							write(*,*)'Angulo theta<10graus ',ilat,jlat
! 							return
						elseif(u .le. xl_el .or. v .le. xl_el)then
							viavel = .false.
							write(*,*)'PC mto proximos',ilat,jlat
							return
						endif
					endif

				enddo
		    
				jlatacum=jlatacum + NUMPCI(jinc)     
			enddo

		enddo

		ilatacum = ilatacum + NUMPCI(inc)
	enddo
	    
	return 
	end
	

	function theta(Ax,Ay,Bx,By,Cx,Cy,u,v)

	double precision Ax,Ay,Bx,By,Cx,Cy
	double precision theta,u1,u2,v1,v2,u,v

	u1 = Bx-Ax
	u2 = By-Ay
	v1 = Cx-Bx
	v2 = Cy-By
	u = sqrt(u1**2 + u2**2)
	v = sqrt(v1**2 + v2**2)


	theta = acos(abs(u1*v1+u2*v2)/(u*v))
	
! 	write(*,*)'u1,u2,u',u1,u2,u,'v1,v2,v',v1,v2,v,'theta ',
!      &			acos((u1*v1+u2*v2)/(u*v))


	return
	end

	
	function intersect(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,u,v)
	
	logical intersect,collinear,Xor,left
	double precision Ax,Ay,Bx,By,Cx,Cy,Dx,Dy
	double precision u1,u2,v1,v2,u,v
	
	if (collinear(Ax,Ay,Bx,By,Cx,Cy) .or. 
     &	    collinear(Ax,Ay,Bx,By,Dx,Dy) .or. 
     &      collinear(Cx,Cy,Dx,Dy,Ax,Ay) .or. 
     &	    collinear(Cx,Cy,Dx,Dy,Bx,By))then
	    
	    intersect = .false.
	else
	    intersect= (Xor( left(Ax,Ay,Bx,By,Cx,Cy), 
     &                       left (Ax,Ay,Bx,By,Dx,Dy) ) .and.  
     &			Xor( left(Cx,Cy,Dx,Dy,Ax,Ay), 
     &			     left(Cx,Cy,Dx,Dy,Bx,By) ) )
	endif

	u1 = Bx-Ax
	u2 = By-Ay
	v1 = Dx-Cx
	v2 = Dy-Cy
	u = sqrt(u1**2 + u2**2)
	v = sqrt(v1**2 + v2**2)

	
	return
	end
	
	
	function collinear(Ax,Ay,Bx,By,Cx,Cy)

	logical collinear
	double precision Ax,Ay,Bx,By,Cx,Cy,area
		
	if( area(Ax,Ay,Bx,By,Cx,Cy) .eq. 0.)then
	    collinear = .true.
	else
	    collinear = .false.
	endif

	return
	end

	function area(Ax,Ay,Bx,By,Cx,Cy)

	
	double precision Ax,Ay,Bx,By,Cx,Cy
	double precision area
	double precision area_parc

	area =  area_parc(Ax,Ay,Bx,By) +
     &		area_parc(Bx,By,Cx,Cy) +
     &		area_parc(Cx,Cy,Ax,Ay)



! 	area = ( ax * by - ay * bx +
!      &           ay * cx - ax * cy +
!      &           bx * cy - cx * by ) !retorna area*2
     
!      	write(*,*)'areas:', area1,area2
	return
	end
	
	function left(Ax,Ay,Bx,By,Cx,Cy)
	
	double precision Ax,Ay,Bx,By,Cx,Cy,area
	logical left
	
	if(area(Ax,Ay,Bx,By,Cx,Cy) .gt. 0.)then
	    left = .true.
	else
	    left = .false.
	endif
     
	return
	end
	
! 	logical function Xor(x,y)
! 	
! 	logical x,y
! 	
! 	if((x .eqv. .false. .and. y .eqv. .true.) .or. 
!      &	   (x .eqv. .true. .and. y .eqv. .false.))then
! 	    Xor = .true.
! 	else
! 	    Xor = .false.
! 	endif
! 	
! 	return
! 	end
	    
	

	subroutine calcDuDQ(dudQ,DB,PAR,numpc,numpar,N_eq,NFIX,NNP)
	
	include 'param_dim.inc'
	double precision dudQ(N_eq,numpar),DB(N_eq,NNP*2),PAR(numpar)
	integer PCnnp(nparg_max) ! conta qtos nohs parametros cada ponto de controle influencia
	integer PCnp(nparg_max,NX) ! guarda quais nohs parametros cada ponto de controle influencia
	double precision PCgama(nparg_max,NX) ! guarda os pesos com que cada coord dos PC influencia cada noh parametro
	double precision PCgamaSx(nparg_max,NX),PCgamaSy(nparg_max,NX)! guarda os pesos com que cada parametro S dos PCs influencia cada Coordenada X ou Y dos nohs parametros
	integer cjNPx,cjNPy,jNP
	
	common/dadossensib/PCnnp,PCnp,PCgama,PCgamaSx,PCgamaSy

        do j=1,numpar
	   do i=1,N_eq
		dudQ(i,j) = 0.
	    enddo
	enddo


	do k=1,numpc !ponto de controle -----Q eh cada coordenada dos PC

            do j=1,PCnnp(k) !loop sobre os nohs que o PC k influencia
                jNP = PCnp(k,j)
                cjNPx = (jNP - NFIX)*2 -1 !coluna de DB referente ao x do noh parametro JNP
                cjNPy = (jNP - NFIX)*2    !coluna de DB referente ao y do noh parametro JNP

	        do i=1,N_eq !loop sobre as derivadas dos valores incognitos

		    if(numpar .eq. 3*numpc)then
			dudQ(i,k*3-2) = dudQ(i,k*3-2) + DB(i,cjNPx)*PCgama(k,j)  !dudQ de Qx
			dudQ(i,k*3-1) = dudQ(i,k*3-1) + DB(i,cjNPy)*PCgama(k,j)  !dudQ de Qy
			dudQ(i,k*3  ) = dudQ(i,k*3  ) + (DB(i,cjNPx)*PCgamaSx(k,j)+ 
     &				DB(i,cjNPy)*PCgamaSy(k,j))/
     &				(3.1416*(PAR(k*3)**2+1.0)) !dudQ de Qs
		    else if(numpar .eq. 2*numpc)then
			dudQ(i,k*2-1) = dudQ(i,k*2-1) + DB(i,cjNPx)*PCgama(k,j)  !dudQ de Qx
			dudQ(i,k*2  ) = dudQ(i,k*2  ) + DB(i,cjNPy)*PCgama(k,j)  !dudQ de Qy
		    else if(numpar .eq. numpc)then
			dudQ(i,k    ) = dudQ(i,k    ) + (DB(i,cjNPx)*PCgamaSx(k,j)+ 
     &				DB(i,cjNPy)*PCgamaSy(k,j))/
     &				(3.1416*(PAR( k )**2+1.0)) !dudQ de Qs
		    endif
		enddo
	    enddo
	enddo


	
	return
	end
	

















*************************************************************************
*                                                                       *
*  Rotina para solucao de sistemas de equacoes lineares nao-simetricos  *
*                                                                       *
*  pelo metodo iterativo Generalized Minimum Residual - GMRES.          *
*                                                                       *
*  Esta implementacao se baseia no artigo de SHAKIB,F. , HUGES,J.R. ,   *
*                                                                       *
*  JOHAN,Z."A Multi-Element Group Preconditioned GMRES Algorithm for    *
*                                                                       *
*  nonsymetric Systems Arrising in Element Analysis" , Computer Methods *
*                                                                       *
*  in Applied Mechanics and Engineering, No 75 ,pags. 415-456, North    *
*                                                                       *
*  Holland, 1989.                                                       *
*                                                                       *
*************************************************************************



        SUBROUTINE GMRES(A,B,XIS,XISA,U,EBAR,EEBAR,HBAR,YIP,CE,
     &                 ES,ETOL,N,NB,K,LMAX,ICONV,ITER,IRP,IWR)


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)



      DIMENSION A(N,N),B(N)
      DIMENSION XIS(N),XISA(N),U(N,NB+1),HBAR(NB+1,NB+1)
      DIMENSION YIP(NB+1),EBAR(NB+1),EEBAR(NB+1),CE(NB+1),ES(NB+1)
        
        dimension baux(1000)

C      COMMON /A/ D(2,2),XI(6,3),W(6,3),IDUP(500),INC(500,2),C(500),
C     *           S(500,3),ISYM(500),X(500),Y(500),IFIP(1000),    
C     *           A(1000,1000),P(1000),B(1000)C
C
C
C      COMMON /GMRS/ XIS(1000),U(1000,91),EBAR(91),HBAR(91,91),
C     *          YIP(91),CE(91),ES(91),           XISA(1000)


C...Descricao das Variaveis

C   N      => Numero de incognitas do sistema.
C   A(N,N) => Matriz de coeficientes do sistema.
C   B(N)   => Vetor independente do sistema.
C   ETOL   => Coeficiente arbitrado para estipular a tolerancia.
C   LMAX   => Numero maximo de iteracoes do GMRES.
C   EPSON  => Tolerancia.
C   XIS(N)   => Vetor que armazena a solucao aproximada.
C   U(N,K) => Matriz que armazena a base ortogonalizada para Z.
C   EBAR(K+1) => Vetor que vai armazenar o residuo na posicao k+1.
C   HBAR(K+1,K+1) => Matriz do sistema que resolvido vai gerar os coef.
C                    das bases ortogonalizadas.
C   YIP(K) => Vetor que acumula os coeficientes das bases ortogonais.




            r=0.
            do i=1,N
                    r = r + b(i)**2
            enddo
            write(*,*)'Antes solve gmres |b|= ',r


C...Inicializacao

       K1= K+1

C......Calcula a norma de B e zera X.
       XNORMX= 0.D0
       XNORMB= 0.D0
       DO 10 IAUX=1,N
          XNORMB= XNORMB + B(IAUX)*B(IAUX)
          XIS(IAUX)= 0.D0
          XISa(IAUX)= 0.D0
 10    CONTINUE
       XNORMB= DSQRT( XNORMB )

       EPSON= ETOL*XNORMB


c       WRITE (6,1000) EPSON
 1000  FORMAT(/,10X,'Norma do residuo admissivel = ',E12.6,/)

c       write(6,1010) XNORMB
 1010  format(/,5x,'Norma inicial do residuo = ',E12.6,/)


C...Ciclos do GMRES

      DO 20 L=1,LMAX
         ll=l

c         !  write(*,1030)
         IF(IRP.GE.1) THEN
           write(IWR,1030) 
 1030      format(/,10X,'Residuo nas iteracoes:'/
     *            /,17x,'Relativo',9x,'Absoluto',12x,
     *            '!!x!!',14x,'dr/dx',/)
         END IF

C...Calculo do primeiro vetor da base ortogonalizada de Z.

C......Zera U(N,K+1)
       DO 25 JAUX=1,K1
          DO 25 IAUX=1,N
             U(IAUX,JAUX)= 0.D0
 25     CONTINUE



C.......Calcula u1= b - A.X

        DO 30 JAUX=1,N
          DO 30 IAUX=1,N

              U(IAUX,1)= U(IAUX,1) - A(IAUX,JAUX)*XIS(JAUX)
 30     CONTINUE

        DO 35 IAUX=1,N

          U(IAUX,1)= U(IAUX,1) + B(IAUX)
 35     CONTINUE

C.......Calcula a norma de U e acumula em EBAR(1).

        EBAR(1)= 0.D0
        DO 50 IAUX=1,N
          EBAR(1)= EBAR(1) + U(IAUX,1)*U(IAUX,1)
 50     CONTINUE
        EBAR(1)= DSQRT( EBAR(1) )

        DO 60 IAUX=1,N
           U(IAUX,1)= U(IAUX,1) / EBAR(1)
 60     CONTINUE


C...Iteracao do GMRES.

        DO 70 I=1,K

           II=I
           I1=I+1

C.........Calcula U(i+1)= A.U(i)

          DO 80 JAUX=1,N
            DO 90 IAUX=1,N

              U(IAUX,I1)= U(IAUX,I1) + A(IAUX,JAUX)*U(JAUX,I)
 90         CONTINUE
 80       CONTINUE



C...Iteracao do algoritmo de ortogonalizacao de Gram-Schmidt modificado.

          DO 100 J=1,I

            HBAR(J,I)=0.D0
            DO 110 IAUX=1,N
               HBAR(J,I)= HBAR(J,I) + U(IAUX,I1)*U(IAUX,J)
 110        CONTINUE

            DO 120 IAUX=1,N
               U(IAUX,I1)= U(IAUX,I1) - HBAR(J,I)*U(IAUX,J)
 120        CONTINUE

 100      CONTINUE


C.........Calcula a norma de U(i+1).

          HBAR(I1,I)=0.D0
          DO 130 IAUX=1,N
            HBAR(I1,I)= HBAR(I1,I) + U(IAUX,I1)*U(IAUX,I1)
 130      CONTINUE

          HBAR(I1,I)= DSQRT( HBAR(I1,I) )

C.........Normaliza o vetor U(i+1).

          DO 140 IAUX=1,N
             U(IAUX,I1)= U(IAUX,I1)/HBAR(I1,I)
 140      CONTINUE


C.........Algoritmo Q-R

          IF(I.GT.1) THEN
            DO 150 J=1,I-1
              HBJI=   CE(J)*HBAR(J,I) + ES(J)*HBAR(J+1,I)
              HBJ1I= -ES(J)*HBAR(J,I) + CE(J)*HBAR(J+1,I)
              HBAR(J,I)= HBJI
              HBAR(J+1,I)= HBJ1I
 150        CONTINUE
          END IF


          R= HBAR(I,I)*HBAR(I,I) + HBAR(I1,I)*HBAR(I1,I)
          R= DSQRT(R)

          CE(I)= HBAR(I,I)/R
          ES(I)= HBAR(I1,I)/R

          HBAR(I,I)= R

          HBAR(I1,I)= 0.D0

          EBAR(I1)= -ES(I)*EBAR(I)
          EBAR(I)=   CE(I)*EBAR(I)

C........ Fim do algoritmo Q-R.

C...Teste de convergencia
          XNORMR= DABS(EBAR(I1))
          XNORXA= XNORMX
          inda=0

          if(i.eq.k) inda=1
          CALL CALCX(XIS,XISA,U,YIP,HBAR,EBAR,EEBAR,inda,I,N,NB,XNORMX)

C          CALL CALCX(inda,I,N,XNORMX)
           DRDX= (XNORMX-XNORXA)/(XNORMR-XNORRA)
          iter= (L-1)*K+I
          IF (ITER.NE.1) THEN
c            !  write(*,1060) iter,XNORMR/xnormb,XNORMR,XNORMX,DRDX
 1060       format(5x,i3,6X,e12.6,5X,e12.6,5X,E16.10,5X,E12.6)
            IF (IRP.GE.1) THEN
              write(IWR,1060) iter,XNORMR/xnormb,XNORMR,XNORMX,DRDX
            END IF
          ELSE
c            !  write(*,1061) iter,XNORMR/xnormb,XNORMR,XNORMX
 1061       format(5x,i3,6X,e12.6,5X,e12.6,5X,E16.10)
            IF (IRP.GE.1) THEN
              write(IWR,1061) iter,XNORMR/xnormb,XNORMR,XNORMX
            END IF
          END IF

          XNORRA= DABS(EBAR(I1))
          IF(DABS(EBAR(I1)).LE.EPSON) GOTO 1

 70     CONTINUE


C...Calculo dos Y por retro substituicao.

 1      I= II
c        CALL CALCX(I,N,XNORMX)
C...Teste de convergencia.

        IF(DABS(EBAR(I1)).LE.EPSON) then
          ICONV =1
          GOTO 2
        END IF

 20   CONTINUE

 2    L=LL

        if(iconv.eq.1)then
          iter= (L-1)*K+I

c          !  write(*,1070)I,L,iter
c          write(IWR,1070)I,L,iter
 1070     format(//,5x,'Convergencia atingida no passo:',i3,/,
     *              5X,'                      do ciclo:',i3,/,
     *              5X,'  ( ',I4,' iteracoes )',/)

        else
          !  write(*,1090)l
          write(IWR,1090)l
 1090     format(//,5x,'NAO foi atingida a convergencia em ',
     *              i3,' ciclos.')
        end if

C.......Calcula u1= b - A.X
        DO 195 IAUX=1,N
           U(IAUX,1)= 0.D0
 195    CONTINUE

        DO 200 JAUX=1,N
          DO 200 IAUX=1,N
              U(IAUX,1)= U(IAUX,1) - A(IAUX,JAUX)*XIS(JAUX)
 200     CONTINUE

        DO 210 IAUX=1,N
          U(IAUX,1)= U(IAUX,1) + B(IAUX)
 210     CONTINUE


C.......Calcula a norma de U e acumula em EBAR(1).

        RESID= 0.D0
        DO 220 IAUX=1,N
          RESID= RESID + U(IAUX,1)*U(IAUX,1)
 220    CONTINUE
        RESID= DSQRT( RESID )

        WRITE(IWR,1100) RESID
 1100   FORMAT(/,5X,'Norma do residuo calculado por  RESID = b - A.x',
     *        //,5x,'    Norma: ',e12.6,/)


C...ARMAZENA OS RESULTADOS NO VETOR B
        DO 190 IAUX=1,N
          B(IAUX)= XIS(IAUX)
 190    CONTINUE

        do i=1,N
            baux(i) = 0.
        enddo

        do i=1,N
            do j=1,N
                baux(i) = baux(i) + b(j)*A(i,j)
            enddo
        enddo


            r=0.
            do i=1,N
                    r = r + baux(i)**2
            enddo
            write(*,*)'depois solve gmres |A.x|= ',r


            r=0.
            do i=1,N
                    r = r + b(i)**2
            enddo
            write(*,*)'depois solve gmres |b|= ',r

            write(*,*)'iteracoes: ',iter,'equacoes: ',N

!         do i=1,N
!             write(*,*)baux(i)
!         enddo
        write(*,*)

      RETURN
      END


      SUBROUTINE CALCX(XIS,XISA,U,YIP,HBAR,EBAR,EEBAR,inda,
     &                 I,N,NB,XNORMX)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DIMENSION XIS(N),XISA(N),U(N,NB+1),HBAR(NB+1,NB+1)
      DIMENSION YIP(NB+1),EBAR(NB+1),EEBAR(NB+1)


c      COMMON /GMRS/ XIS(1000),U(1000,91),EBAR(91),HBAR(91,91),
c     *          YIP(91),CE(91),ES(91),        xisa(1000)

C      dimension eebar(NB+1)

        I1= I+1
        IM1= I-1

        do 10 iaux=1,nb+1
            eebar(iaux)= ebar(iaux)
 10     continue

        YIP(I)= eEBAR(I)/HBAR(I,I)

        DO 160 IAUX=1,IM1

          KAUX= I-IAUX
          KAUX1= KAUX+1

          DO 170 JAUX=KAUX1,I
            eEBAR(KAUX)= eEBAR(KAUX) - HBAR(KAUX,JAUX)*YIP(JAUX)
 170      CONTINUE
          YIP(KAUX)= eEBAR(KAUX)/HBAR(KAUX,KAUX)
 160    CONTINUE

        do 175 iaux=1,n
           xis(iaux)=xisa(iaux)
 175    continue

        DO 180 JAUX=1,I
          DO 180 IAUX=1,N
             XIS(IAUX)= XIS(IAUX) + YIP(JAUX)*U(IAUX,JAUX)
 180      CONTINUE

C.......CALCULA A NORMA DA SOLUCAO.
        XNORMX=0.D0
        DO 190 IAUX=1,N
           XNORMX= XNORMX + XIS(IAUX)*XIS(IAUX)
 190    CONTINUE
        XNORMX= DSQRT(XNORMX)

C.......Atualiza a solucao.
       if (inda.ne.0) then 
        do 200 iaux=1,n
           xisa(iaux)=xis(iaux)
 200    continue
        end if

        RETURN
        END



      	SUBROUTINE sensib_OUTPT_OT(dudQ,dudQOT,numpar,numpc,ldfjac,
     &                    I_ELET,IDCASO,IC,NEPE,NCasos,N_eq)

	INCLUDE 'param_dim.inc'		
		
	double precision dudQ(N_eq,numpar),dudQOT(ldfjac,numpar)
	double precision aux(npar_max)
	
	INTEGER IDCASO(NCasosX,3),I_ELET(16,NELEX)

	IELET1= IDCASO(IC,1)
	IELET2= IDCASO(IC,2)
	IELET3= IDCASO(IC,3)

! 	IR0= (IC-1)*(16-3)
	IR0= (IC-1)*(16-2)
	IRES= 0
! 	write(*,*)'dudqOT'	
	DO IE= 1,16
		IF( (IE.NE.IELET1).AND.
!      &       	    (IE.NE.IELET3).AND.  ! A partir da modif da funcao objetivo usa-se o potencial do elet de ref
     &              (IE.NE.IELET2) ) THEN
			IRES = IRES + 1

			INOF = I_ELET(IE,1+NEPE/2)

			if(IE.NE.IELET3)then
				do j=1,numpar
					dudQOT(IR0+IRES,j) = dudQ(INOF,j)  !derivadas dos valores incognitos
				enddo
			else
				do j=1,numpar
					dudQOT(IR0+IRES,j) = 0. !derivadas do valor prescrito
				enddo		
			endif

		ENDIF
	ENDDO



!  Consequencia da regra da cadeia J(i,j) = DuDQOT(i,j) + 1/14* sum_k (DuDQOT(k,j)),k=1,14
! Dc/Dz_i = 1/14 para qq z_i



	do j = 1,numpar
		aux(j) = 0.
		do k = (IR0+1),(IR0+16-2)
			aux(j) = aux(j) + dudQOT(k,j)
		enddo
	enddo

	do i=(IR0+1),(IR0+16-2)
		do j=1,numpar
			dudQOT(i,j) = dudQOT(i,j) - (1./14.)*aux(j)
		enddo
	enddo

! 	do i=IR0+1,IR0+16-2
! 		write(4,2)(dudQOT(i,j),j=1,numpar)
! 	enddo




 1	format(i10,15f10.5)
 2	format(8f10.5)
	return
	end

      SUBROUTINE OUTPT_OT(FI,DFI,RCONT,RCEX,RESID,I_ELET,
     &                  IDCASO,IC,NEPE,NCasos,FDEX,FROT)


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      INCLUDE 'param_dim.inc'

      DIMENSION FI(NX),DFI(NX)
      DIMENSION RCONT(NRCX),RCEX(NRCX)
	  real*8 c

      INTEGER IDCASO(NCasosX,3),I_ELET(16,NELEX)
      INTEGER FDEX,FROT
      real resid

      IELET1= IDCASO(IC,1)
      IELET2= IDCASO(IC,2)
      IELET3= IDCASO(IC,3)

!       IR0= (IC-1)*(16-3)
      IR0= (IC-1)*(16-2)
      IRES= 0

      DO IE= 1,16
         IF( (IE.NE.IELET1).AND.
!      &       (IE.NE.IELET3).AND.   A partir da modif da funcao objetivo usa-se o potencial do elet de ref
     &       (IE.NE.IELET2) ) THEN
            IRES= IRES + 1
            INOF= I_ELET(IE,1+NEPE/2)
	    RCONT(IR0+IRES)= FI(INOF)
         ENDIF
      ENDDO



! ! ! Alteracao do residuo por uma constante "c"
	c = 0.
	do i=IR0+1,IR0+IRES
		c = c + (RCONT(i) - RCEX(i))
	enddo

	c = c / IRES !calculo da constante "c"
! 	write(*,*) 'c: ',c,ires
	
	do i=IR0+1,IR0+IRES !alteracao das diferencas de potencial
		RCONT(i) = RCONT(i) - c
	enddo
! ! ! 


      IF(IC.EQ.NCASOS) THEN

!          NRC= NCASOS*(16-3)
         NRC= NCASOS*(16-2)

         RESID= 0.D0
         DO I=1,NRC
            RESID= RESID + (RCONT(I)-RCEX(I))*(RCONT(I)-RCEX(I))
         ENDDO
         RESID= 0.5 * RESID

         WRITE(FROT,1010) (RCEX(I),I=1,NRC)

         WRITE(FROT,1010) ((RCONT(I)),I=1,NRC)
 1010    FORMAT(8(14E12.4,/))

         WRITE(FROT,1020) RESID
 1020    FORMAT(5X,'RESIDUO :',E15.5)
      
      ENDIF


     
      RETURN
      END


	subroutine input_opt(ftol,xtol,gtol,maxfev,numpar,
     &                 NRC,metodo,mode,nprint,factor,epsfcn)
        include 'param_dim.inc'	
	COMMON /cFDAT/ FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES
	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA	
	CHARACTER*15 FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES
	integer FOPT,FDAT
	
	double precision p(npar_max),s(npar_max)
	double precision ftol,xtol,gtol
	integer maxfev,numcc,numpar,NRC,mode,metodo,nprint
	double precision factor,epsfcn
	double precision PGLIM(nparg_max,3),CINC(ninc_max),C0
	integer VARIA(ninc_max)	
	integer NNINC(ninc_max),NUMPI(ninc_max),NUMPCI(ninc_max)
	common/geom/p,s,numpc

	FILE_DAT='bacia_ext00.dat'
	
	WRITE(*,*) 'Arquivo cont.ext.: bacia_ext00.dat ? (0=>NAO)'
	READ(*,*) ISIM
	
	IF(ISIM.EQ.0) THEN
            WRITE(*,'(a/)')' Entre arq.DADOS GEOM.FIXOS(10 char.bas.):'
            READ(*,'(A)') FILE_DAT
	END IF
	
	WRITE(*,'(a/)')' Arq.. MEDIDAS NOS ELETRODOS (10 char. bas.):'
	READ(*,'(A)') FILE_DEX
	
	write(*,*) 'Entre com ARQUIVO DE DADOS DE OTIMIZACAO:'
	READ(*,*) FILE_OPT


	OPEN(4,file=file_opt//'.out',status='new')

!  SAIDA 
	WRITE(4,*)'Arq.DADOS GEOM.FIXOS(10 char.bas.):'
	WRITE(4,*) FILE_DAT

	WRITE(4,*)'Arq.. MEDIDAS NOS ELETRODOS (10 char. bas.):'
	WRITE(4,*) FILE_DEX
	
	WRITE(4,*) 'ARQUIVO DE DADOS DE OTIMIZACAO:'
	WRITE(4,*) FILE_OPT
! 
	
	FOPT= 11
	OPEN (FOPT,FILE=FILE_OPT,STATUS='OLD')


	READ(FOPT,*) ftol,xtol,gtol,maxfev

	write(4,*)
	write(4,*)'   Parametros do alg. de otimizacao'
	write(4,*)
	write(4,*)'   ftol: ',ftol
	write(4,*)'   xtol: ',xtol
	write(4,*)'   gtol: ',gtol
	write(4,*)'   maxfev: ',maxfev

	READ(FOPT,*) metodo
	if(metodo .eq. 1)then
		READ(FOPT,*) mode,factor,nprint,epsfcn
		write(4,*)'   mode: ',mode
		write(4,*)'   factor: ',factor
		write(4,*)'   nprint: ',nprint
		write(4,*)'   epsfcn: ',epsfcn
	else if(metodo .eq. 2)then
		READ(FOPT,*) mode,factor,nprint
		write(4,*)'   mode: ',mode
		write(4,*)'   factor: ',factor
		write(4,*)'   nprint: ',nprint
	endif

	READ(FOPT,*) NINCL


	write(4,1020) NINCL
 1020   format( //,3x,'Numero de inclusoes: ',i5,/)

	NUMPAR = 0
	NUMPC = 0
	do i=1,nincl
		read(fopt,*)NUMPCI(i),NUMPI(i),cinc(i),VARIA(i)
		write(4,*)'   N.PCs: ',NUMPCI(i),'  N.Params: ',NUMPI(i),
     &			'  Condut.: ',cinc(i),'  Varia?: ',VARIA(i)
		NUMPAR = NUMPAR + NUMPI(i)
		NUMPC = NUMPC + NUMPCI(i)
	enddo

	if(NUMPAR .NE. NUMPC .and. 
     &		NUMPAR .NE. NUMPC*2 .and. 
     &		NUMPAR .NE. NUMPC*3)then
		write(*,*) 'ERRO: num. de param. ou pontos de controle'
		write(4,*) 'ERRO: num. de param. ou pontos de controle'
		stop
	endif

	write(4,*) 
	if(NUMPAR .EQ. NUMPC)then
		write(4,*) 'variaveis de otimizacao: parametro S'
	else if(NUMPAR .EQ. 2* NUMPC)then
		write(4,*) 'variaveis de otimizacao: coordenadas'
	else if(NUMPAR .EQ. 3*NUMPC)then
		write(4,*) 'variaveis de otimizacao: coord. e param. S'
	endif
		write(4,*)

	i=0
	do iinc = 1, nincl
		do ipi= 1, NUMPCI(iinc)
			i=i+2
			READ(FOPT,*) P(i-1),P(i),S(i/2)
		enddo
	enddo

	WRITE(4,*) ' Pontos de Controle Iniciais'
	WRITE(4,*) ' inclusao       X          Y'

	i=0
	do iinc = 1, nincl
		do ipi= 1, NUMPCI(iinc)
			i=i+2
			write(4,1000)iinc,P(i-1),P(i),S(i/2)	
		enddo
	enddo
 1000   FORMAT(I5,3F15.4)


	FDAT = 22
	open(FDAT,file = FILE_DAT)
	READ(FDAT,150)TITULO
  150 	FORMAT(18A4)
	READ(FDAT,*)NFIX,NEPE,NCasos,ITCC, VCC
! 	NRC = Ncasos*(16-3)
	NRC = Ncasos*(16-2)
	close(FDAT)
	
	return
	end




      SUBROUTINE INPUT0(X,Y,RCEX,VCC,INC,NFIX,NEPE,I_ELET,NCasos,
     &                  IDCASO,ITCC,NG,NE)
                        

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      INCLUDE 'param_dim.inc'

      INTEGER       FDAT,FDCC,FDEX,FRES,FROT
      COMMON /narq/ FDAT,FDCC,FDEX,FRES,FROT



      COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA

      

      character TITULO*72
      DIMENSION X(NX),Y(NX),INC(NX,2)
      DIMENSION RCEX(NRCX)

      DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max),NUMPCI(ninc_max)
      integer VARIA(ninc_max)
      DIMENSION PGLIM(nparg_max,3)      

      INTEGER IDCASO(NCasosX,3),I_ELET(16,NELEX)

C
C N=   Numero de elementos de contorno.
C L=   Numero de pontos internos onde a funcao e calculada.
C
      WRITE(FRES,100)
  100 FORMAT(' ',80('*'))
C
C Leitura do Titulo do problema.
C
      READ(FDAT,150) TITULO
  150 FORMAT(A72)
      WRITE(FRES,250) TITULO
!        write(*,250) titulo

  250 FORMAT(5X,A72)
C
C Leitura dos parametros basicos.
C
      READ(FDAT,*)NFIX, NEPE, NCasos, ITCC, VCC, C0
 200  FORMAT(4I10,2F10.4)

C     Inicialmente NEPE deve ser considerado impar e o 
C     potencial do eletrodo a ser considerado ser� o do
C     elemento do meio.

C     Numero de elementos fixos igual ao nimero de nos 
      NEFIX= NFIX

C Numero de pontos de integra��o
      NG=8

      WRITE(FRES,300) NFIX,NFIX,NG,NEPE,NCasos,ITCC,VCC,C0
  300 FORMAT(//2X,'DADOS'//,
     *  5x,'NUMERO DE NOS PARA DEFINICAO DA GEOMETRIA =',I4,/,
     *  5X,'NUMERO DE ELEMENTOS DE CONTORNO           =',I4,/,
     *  5X,'NUMERO DE PONTOS DE GAUSS                 =',I4,/,
     *  5X,'NUMERO DE ELEMENTOS POR ELETRODO          =',I4,/,
     *  5X,'NUMERO DE CASOS DE CARREGAMENTO           =',I4,///,
     *  5X,'TIPO DE SOLICITACAO (1=> Corrente)        =',I4,/,
     *  5X,'VALOR DA SOLICITACAO                      =',F12.4,/,
     &  5x,'CONDUTIVIDADE DO MEIO                     =',F12.4,/)

      IF (NFIX.GT.NX) STOP 'NUM. INCOGNITAS MAIOR QUE O DIMENSIONADO' 

C
C Leitura das coordenadas dos pontos extremos dos elementos de
C contorno nos vetores X e Y.
C
      WRITE(FRES,500)
  500 FORMAT(//2X,'COORDENADAS DOS PONTOS EXTREMOS DOS ELEMENTOS DE '
     *       ,'CONTORNO FIXOS',//4X,'PONTOS',10X,'X',18X,'Y')
      DO 10 I=1,NFIX
         READ(FDAT,*)k,X(k),Y(k)
 10   CONTINUE

      DO I=1,NFIX
         WRITE(FRES,700)I,X(I),Y(I)
  700    FORMAT(5X,I3,2(5X,E14.7))
	END DO

C
C Leitura das incidencias dos elementos.
C
      DO I=1,NFIX
         READ(FDAT,*) K,INC(K,1),INC(K,2)
C 750     FORMAT(*)
      END DO

      WRITE(FRES,760)
 760  FORMAT(//2X,'INCIDENCIA DOS ELEMENTOS FIXOS',//4X,'ELEMENTO',5X,
     &            'NO INICIAL',5X,'NO FINAL',5x,'COMPRIMENTO')

      DO I=1,NFIX
         dx= x(inc(i,2)) - x(inc(i,1))
         dy= y(inc(i,2)) - y(inc(i,1))
         xl = sqrt(dx*dx + dy*dy)
         WRITE(FRES,770) I, INC(I,1), INC(I,2), xl
 770     FORMAT(6X,I3,12X,I3,9X,I3,5X,F12.4)
      END DO



C
C Leitura das condicoes de contorno.
C FI(I)= Valor do potencial no noh se KODE=0 .

      
C     Le E ESCREVE os nos que compoem os eletrodos

      DO I = 1,16
         READ(FDCC,*) (I_ELET(I,J),J=1,NEPE)
      ENDDO


      WRITE(FRES,780) 
 780  FORMAT(//,2X,'DISCRETIZACAO DOS ELETRODOS',//,
     &          2X,'ELETRODO',5X,'ELEMENTOS',/)     

      DO I=1,16
         WRITE(FRES,790) I,(I_ELET(I,J),J=1,NEPE)
 790     FORMAT(2X,I5,10X,3I5)
      ENDDO


C     Le a definicao dos casos de carga em funcao dos eletrodos
C     eletrodo de entrada
C     eletrodo de saida 
C     eletrodo de referencia (potencial 0)

      DO I =1,NCasos
         READ(FDCC,*) (IDCASO(i,j),j=1,3)
      ENDDO

C     Impressao dos Casos de Carga

      WRITE(FRES,800)
 800  FORMAT(/,'  Definicao dos Casos de Carga ',//,
     &         '  Caso     Eletrodo    Eletrodo   Eletrodo',/,
     &         '           Entrada      Saida     Referencia',/)
      DO I =1,NCasos
         WRITE(FRES,1020) I, (IDCASO(i,j),j=1,3)
 1020    FORMAT( I5,3i12)
      ENDDO

!       NRC= NCASOS*(16-3)
      NRC= NCASOS*(16-2)

      READ(FDEX,*)(RCEX(I),I=1,NRC) 
      
      read(fdex,*) !linha vazia
      read(fdex,*)snr
      open(1,file='snr.tmp')
      write(1,*)snr
      close(1)


      WRITE(FRES,1030) (RCEX(I),I=1,NRC)
 1030 FORMAT(21(5E12.4,/))

      RETURN
      END


	SUBROUTINE GERACONT(PARG,XL_EL,X,Y,NUMPC,FRES,INC,NFIX,
     &			NN,N_eq,npa0) 

	implicit double precision(a-h,o-z)
	
	INCLUDE 'param_dim.inc'
	
	common /NPAR/ NPAR
	
	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA
! 	common/plotcav/NNINCaux,NINCLaux
	common /triang/ct
	real ct(nparg_max/2,6)
	DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max),NUMPCI(ninc_max)
	integer VARIA(ninc_max)
	
	DIMENSION NNINCaux(ninc_max)
	dimension xlac(ntmx-1),t(ntmx),xt(ntmx),yt(ntmx),
     &            xk(NPC_MAX+3),yk(NPC_MAX+3),sk(NPC_MAX+3),  
     &            tel(NX),sel(NX),
     &            X(NX),Y(NX),INC(NX,2)

	real*8 PARG(nparg_max)
	integer fres


	integer PCnnp(nparg_max) ! conta qtos nohs parametros cada ponto de controle influencia
	integer PCnp(nparg_max,NX) ! guarda quais nohs parametros cada ponto de controle influencia
	double precision PCgama(nparg_max,NX) ! guarda os pesos com que cada PC influencia cada noh parametro
	double precision PCgamaSx(nparg_max,NX),PCgamaSy(nparg_max,NX)! guarda os pesos com que cada parametro S dos PCs influencia cada Coordenada X ou Y dos nohs parametros

	common/dadossensib/PCnnp,PCnp,PCgama,PCgamaSx,PCgamaSy


C     GERA UMA SPLINE PARA CADA INCLUSAO E ACUMULA O NUMERO DE NOS 
C     DE CADA SPLINE REFERENTE A INCLUS�O i NO VETOR NNINC(i)


C     Inicializa Numero de Equa��es e numero de nos acumulado.
	N_Eq = NFIX
	NN = NFIX
	nna = NN
	nna0= NN


	
C     Parametros Geometricos
	npa0= 0
	
	do iincl=1,nincl

C     Numero de parametros geometricos da inclusao e 
            npi= numpci(iincl)*3	



C     Numero de parametros geometricos 
C     acumulado das inclusoes anteriores.
            npa0= npa0 + npi
    
            !  write(*,*) 'Parametros da Spline:',npi/2
            !  write(*,*) 'Par.       X         Y'
            write(fres,*) 'Parametros da Spline:',npi/3
            write(fres,*) 'Par.       X         Y'
    
            !  write(*,1100) (i,parg(2*i-1), parg(2*i), i=1,npi/2)
            write(fres,1100)(i,PARG(3*i-2),PARG(3*i-1),
     &              PARG(3*i),i=1,npi/3)
 1100       format( i5,3f10.4)

C     Calcula o numero total de pontos para a
C     definicao da XSPLINE e os respectivos t.

            nk=npi/3+1
            
            xk(1)=PARG((npa0-npi)+npi-2)
            yk(1)=PARG((npa0-npi)+npi-1)
            sk(1)=PARG((npa0-npi)+npi  )


            do i= 1,npi/3
                xk(i+1)=PARG((npa0-npi)+i*3-2)
                yk(i+1)=PARG((npa0-npi)+i*3-1)
                sk(i+1)=PARG((npa0-npi)+i*3  )
            end do

            xk(nk+1)=PARG((npa0-npi)+1)
            yk(nk+1)=PARG((npa0-npi)+2)
            sk(nk+1)=PARG((npa0-npi)+3)
    

            xk(nk+2)=PARG((npa0-npi)+4)
            yk(nk+2)=PARG((npa0-npi)+5)
            sk(nk+2)=PARG((npa0-npi)+6)

! 		sk(1)=0
! 		do i=2,nk+1
! 			sk(i)=1
! 		end do
! 		sk(nk+2)=0
          
       
! 		do i=1,nk+2
! 			write(*,*) xk(i),yk(i),sk(i)
! 		end do

c       !  write(*,*) 'nk= ',nk

            nt= ntint*(nk-1) + nk
            dlt= 1.d0/(ntint+1)
            
            
            do it=1,nt
                    t(it)= (it-1)*dlt
            end do
		
            idadossensib = 0 ! nao salva dados uteis para calculo das sensibilidades
            call xsplxy(xk,yk,sk,t,xt,yt,nk,nt,
     &			nna0,nfix,npa0-npi,idadossensib) 


C     Calcula o comprimento dos segmentos da XSPLINE

C     Zera os comprimentos
            do it=1,nt
                    xlac(it)=0.d0
            end do

            do it=1,nt-1
                dltx= (xt(it+1)-xt(it))
                dlty= (yt(it+1)-yt(it))
                xlac(it+1)=xlac(it)+sqrt(dltx*dltx+dlty*dlty)
            end do

!                 write(*,*)nt,xlac(nt),dltx,dlty
		NPAR=int(xlac(nt)/XL_EL)

C     NUMERO DE ELEMENTOS POR INCLUSAO
		NNINC(IINCL)= NPAR
		
		c=1.d0/NPAR
	
		nna= nna0 + NPAR

		N_Eq = N_Eq + 2*npar
c       NE=N

		if(N_Eq.gt.nx) THEN 
			!  write(*,*)'N= ',N_Eq,'NX= ',NX
			stop'N (Neq) maior que NX em geracont '
		ENDIF
		

		do iel=1,(NPAR-1)
			sel(iel)=iel*c
		end do   

C     Calcula os t dos elementos (tel) com base nos SEL
C     (comprimentos relativos prescritos acumulados)

		tel(1)=0.d0
		tel(NPAR+1)=nk-1
		it=1
		do iel=2,NPAR
			xl= sel(iel-1)*xlac(nt)
 1       		if (xl.gt.xlac(it)) then
				it=it+1
				go to 1
			end if
			tel(iel) = dlt*((it-2) + 
     &                 (xl-xlac(it-1))/(xlac(it)-xlac(it-1)))
       		end do

		idadossensib = 1 ! salva dados uteis para calculo das sensibilidades
		call xsplxy(xk,yk,sk,tel,X(nna0+1),Y(nna0+1),nk,NPAR+1,
     & 			nna0,nfix,npa0-npi,idadossensib)	
! 		call xsplxy(xk,yk,sk,tel,X(nna0+1),Y(nna0+1),nk,NPAR+1,
!      & 			nna0,nfix,npa0-npi,idadossensib)	



		do i=1,NPAR
			INC(nna0+i,1)= nna0 + i
			INC(nna0+i,2)= nna0 + i + 1
		end do
         
		INC(nna,2)= nna0 + 1	


      		WRITE(FRES,10)
 10   		FORMAT(//2X,'COORDENADAS DOS PONTOS EXTREMOS DOS ELEMENTOS DE ',
     & 'CONTORNO MOVEIS',//4X,'PONTOS',10X,'X',18X,'Y')

		DO I=nna0+1,nna
			WRITE(FRES,700)I,X(I),Y(I)
 700      		FORMAT(5X,I3,2(5X,E14.7))
		END DO

		WRITE(FRES,20)
 20 		FORMAT(//2X,'INCIDENCIA DOS ELEMENTOS MOVEIS',//4X,'ELEMENTO',5X,
     & 'NO INICIAL',5X,'NO FINAL',5x,'COMPRIMENTO')

		DO I=nna0+1,nna
			dx= x(inc(i,2)) - x(inc(i,1))
			dy= y(inc(i,2)) - y(inc(i,1))
			xl = sqrt(dx*dx + dy*dy)
			WRITE(FRES,710) I, INC(I,1), INC(I,2), xl
 710     		FORMAT(6X,I3,12X,I3,9X,I3,5X,F12.4)
		END DO
  
C     Atualiza o numero de n�s anterior e o numero de equa�oes
		nna0= nna0 + npar

! 		NNINCaux(iincl) = NNINC(iincl) !dado auxiliar para common

c     finaliza o loop sobre as inclusoes
	enddo

! 	do i=1,npa0
! 		write(*,11)(dQdt(i,j),j=1,numpar)
!  11 		format(15f10.3)
! 	enddo
! 	NINCLaux = NINCL !dado auxiliar para common

C     Define o numero final de nos como o numero de nos acumulado
	NN  = nna

	return
	end




      double precision function fpt(p,t)
      implicit double precision(a-h,o-z)
         
          fpt = (10.d0 - p + (2.d0*p-15.d0)*t + (6.d0-p)*t*t )*t*t*t
       return
       end

      double precision function DfptDp(p,t)
      implicit double precision(a-h,o-z)
         
          DfptDp = -t**3 + 2.0*t**4 - t**5

       return
       end

      double precision function DfptDt(p,t)
      implicit double precision(a-h,o-z)
         
          DfptDt = 30.0*t**2 - 3.0*p*t**2 +
     *             4.*(2.*p -15.)*t**3 + 5.0*(6.0 - p)*t**4

       return
       end


	subroutine xsplxy(xk,yk,sk,t,xt,yt,nk,nt,
     &			nna0,nfix,npa0,idadossensib)

	implicit double precision(a-h,o-z)
	INCLUDE 'param_dim.inc'
	dimension xk(nk+3),yk(nk+3),sk(nk+3),t(nt),xt(nt),yt(nt)
	integer PCnnp(nparg_max) ! conta qtos nohs parametros cada ponto de controle influencia
	integer PCnp(nparg_max,NX) ! guarda quais nohs parametros cada ponto de controle influencia
	double precision PCgama(nparg_max,NX) ! guarda os pesos com que cada coordenada dos PCs influencia cada noh parametro
	double precision PCgamaSx(nparg_max,NX),PCgamaSy(nparg_max,NX)! guarda os pesos com que cada parametro S dos PCs influencia cada Coordenada X ou Y dos nohs parametros

	common/dadossensib/PCnnp,PCnp,PCgama,PCgamaSx,PCgamaSy

	if(idadossensib .eq.1)then
		npc = nk - 1	
		do i=npa0/3+1  , npa0/3 + npc
			PCnnp(i)=0 ! zera os contadores de nohs para os PC's desta spline
		enddo

! 		do i=1,nparg_max
! 			do j=1, NX
! 				PCgamaSx(i,j) = 0.
! 				PCgamaSy(i,j) = 0.
! 			enddo
! 		enddo

	endif

	do it=1,nt	

		itk1= 2 + t(it)
		itk2= itk1+1
		
		xk0= xk(itk1-1)  ! cuidado
		xk1= xk(itk1)
		xk2= xk(itk2)
		xk3= xk(itk2+1)  ! cuidado
		
		yk0= yk(itk1-1)  ! cuidado
		yk1= yk(itk1)
		yk2= yk(itk2)
		yk3= yk(itk2+1)  ! cuidado
		
		sk0= sk(itk1-1)  ! cuidado
		sk1= sk(itk1)
		sk2= sk(itk2)
		sk3= sk(itk2+1)  ! cuidado

		Tkp = itk1 + sk1
		Tk1p= itk2 + sk2
		Tk2m= itk1 - sk1
		Tk3m= itk2 - sk2
		
		pkm1= 2.D0*( (itk1-1) - Tkp )*( (itk1-1) - Tkp )
		pk  = 2.D0*(  itk1    - Tk1p)*(  itk1    - Tk1p)
		pkp1= 2.D0*(  itk2    - Tk2m)*(  itk2    - Tk2m)
		pkp2= 2.D0*( (itk2+1) - Tk3m)*( (itk2+1) - Tk3m)

		if(idadossensib .eq.1)then
			Dpkm1DTkp  = -4.0*((itk1-1) - Tkp )
			DpkDTk1p   = -4.0*( itk1    - Tk1p)
			Dpkp1DTk2m = -4.0*( itk2    - Tk2m)
			Dpkp2DTk3m = -4.0*((itk2+1) - Tk3m)
		endif

		t0= (t(it) + 2.d0 - Tkp )/( itk1-1 - Tkp )
		t1= (t(it) + 2.d0 - Tk1p)/( itk1   - Tk1p)
		t2= (t(it) + 2.d0 - Tk2m)/( itk2   - Tk2m)
		t3= (t(it) + 2.d0 - Tk3m)/( itk2+1 - Tk3m)

		if(idadossensib .eq.1)then
			Dt0DTkp  = -1.0 / (  itk1-1 - Tkp  ) + 
     &			(t(it) + 2.d0 - Tkp )/(  itk1-1 - Tkp  )**2
			Dt1DTk1p = -1.0 / (  itk1    - Tk1p) + 
     &			(t(it) + 2.d0 - Tk1p)/(  itk1    - Tk1p)**2
			Dt2DTk2m = -1.0 / (  itk2    - Tk2m) + 
     &			(t(it) + 2.d0 - Tk2m)/(  itk2    - Tk2m)**2
			Dt3DTk3m = -1.0 / ( itk2+1 - Tk3m) + 
     &			(t(it) + 2.d0 - Tk3m)/( itk2+1 - Tk3m)**2
		endif


		if(t(it)+2.d0 .gt. Tkp) then
			A0=0.d0
		else
			A0= fpt(pkm1, t0)
		endif
		A1=  fpt(pk  , t1)
		A2=  fpt(pkp1, t2)

		if(t(it)+2.d0 .lt. Tk3m) then
			A3=0.d0
		else
			A3= fpt(pkp2, t3)
		endif

		if(idadossensib .eq.1)then
			if(t(it)+2.d0 .gt. Tkp) then
				DA0Dsk1=0.d0
			else
				DA0Dsk1= +1.0*(DfptDp(pkm1, t0) * Dpkm1DTkp + 
     &					DfptDt(pkm1, t0) * Dt0DTkp)
			endif

			DA1Dsk2= +1.0*(DfptDp(pk  , t1) * DpkDTk1p + 
     &					DfptDt(pk  , t1) * Dt1DTk1p ) 
			DA2Dsk1= -1.0*(DfptDp(pkp1, t2) * Dpkp1DTk2m + 
     &					DfptDt(pkp1, t2) * Dt2DTk2m)
	
			if(t(it)+2.d0 .lt. Tk3m) then
				DA3Dsk2 =0.d0
			else
				DA3Dsk2 = -1.0*(DfptDp(pkp2, t3) * Dpkp2DTk3m + 
     &					DfptDt(pkp2, t3) * Dt3DTk3m)
			endif
		endif

		xt(it)=(A0*xk0 + A1*xk1 + A2*xk2+ A3*xk3)/(A0+A1+A2+A3)
		yt(it)=(A0*yk0 + A1*yk1 + A2*yk2+ A3*yk3)/(A0+A1+A2+A3)


		if(idadossensib .eq.1 .and. it .lt. nt)then
			! salva dados uteis para o calculo da sensibilidades
			! ip0,ip1,ip2,ip3 sao os 4 pontos de controle envolvidos
			ip3= itk1 + 2 -1 ! o -1 é para relacionar xk com a ordem dos pc dados
			if(ip3 .gt. npc) ip3 = ip3 - npc
			ip2 =ip3 - 1
			if(ip2 .lt. 1) ip2 = ip2 + npc
			ip1 =ip2 - 1
			if(ip1 .lt. 1) ip1 = ip1 + npc	
			ip0 =ip1 - 1
			if(ip0 .lt. 1) ip0 = ip0 + npc	
	! 	 	write(*,*)'it: ',it,'pontos', ip0,ip1,ip2,ip3
	
			ip0 = ip0 + npa0/3 ! passa ip local para ip global
			ip1 = ip1 + npa0/3
			ip2 = ip2 + npa0/3
			ip3 = ip3 + npa0/3
			



		
			PCnnp(ip0)=PCnnp(ip0)+1
			PCnnp(ip1)=PCnnp(ip1)+1
			PCnnp(ip2)=PCnnp(ip2)+1
			PCnnp(ip3)=PCnnp(ip3)+1

			if(npc .eq. 3)then ! neste caso, o mesmo PC foi incrementado duas vezes (no 1o e no 4o xk)
				PCnnp(ip3)=PCnnp(ip3)-1 !desfaz o incremento duplicado
			endif
		
		
			idnoh=nna0+it
! 			if(it .eq. nt)then
! 			    idnoh = nna0+1
! 			endif
			PCnp(ip0,PCnnp(ip0)) = idnoh
			PCnp(ip1,PCnnp(ip1)) = idnoh
			PCnp(ip2,PCnnp(ip2)) = idnoh
			PCnp(ip3,PCnnp(ip3)) = idnoh
		
			PCgama(ip0,PCnnp(ip0)) = A0 / (A0+A1+A2+A3)
			PCgama(ip1,PCnnp(ip1)) = A1 / (A0+A1+A2+A3)
			PCgama(ip2,PCnnp(ip2)) = A2 / (A0+A1+A2+A3)
			PCgama(ip3,PCnnp(ip3)) = A3 / (A0+A1+A2+A3)

			PCgamaSx(ip0,PCnnp(ip0)) = 0.0
			PCgamaSx(ip1,PCnnp(ip1)) = (DA0Dsk1*xk0 + DA2Dsk1*xk2)/
     &			(A0+A1+A2+A3)-(A0*xk0+A1*xk1+A2*xk2+A3*xk3)*
     &			(DA0Dsk1 + DA2Dsk1)/(A0+A1+A2+A3)**2 
			PCgamaSx(ip2,PCnnp(ip2)) = (DA1Dsk2*xk1 + DA3Dsk2*xk3)/
     &			(A0+A1+A2+A3)-(A0*xk0+A1*xk1+A2*xk2+A3*xk3)*
     &			(DA1Dsk2 + DA3Dsk2)/(A0+A1+A2+A3)**2
			PCgamaSx(ip3,PCnnp(ip3)) = 0.0	

			PCgamaSy(ip0,PCnnp(ip0)) = 0.0
			PCgamaSy(ip1,PCnnp(ip1)) = (DA0Dsk1*yk0 + DA2Dsk1*yk2)/
     &			(A0+A1+A2+A3)-(A0*yk0+A1*yk1+A2*yk2+A3*yk3)*
     &			(DA0Dsk1 + DA2Dsk1)/(A0+A1+A2+A3)**2 
			PCgamaSy(ip2,PCnnp(ip2)) = (DA1Dsk2*yk1 + DA3Dsk2*yk3)/
     &			(A0+A1+A2+A3)-(A0*yk0+A1*yk1+A2*yk2+A3*yk3)*
     &			(DA1Dsk2 + DA3Dsk2)/(A0+A1+A2+A3)**2
			PCgamaSy(ip3,PCnnp(ip3)) = 0.0	

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! ! ! ! !  Verificacao por diferencas finitas
! ! ! ! 			deltaS=0.01
! ! ! ! 	
! ! ! ! 			do i=0,3
! ! ! ! 
! ! ! ! 				if(i.eq.0)then
! ! ! ! 					sk0 = sk0 + deltaS
! ! ! ! 					ip = ip0
! ! ! ! 				else if(i.eq.1)then
! ! ! ! 					sk0 = sk0 - deltaS
! ! ! ! 					sk1 = sk1 + deltaS
! ! ! ! 					ip = ip1
! ! ! ! 				else if(i.eq.2)then
! ! ! ! 					sk1 = sk1 - deltaS
! ! ! ! 					sk2 = sk2 + deltaS
! ! ! ! 					ip = ip2
! ! ! ! 				else if(i.eq.3)then
! ! ! ! 					sk2 = sk2 - deltaS
! ! ! ! 					sk3 = sk3 + deltaS
! ! ! ! 					ip = ip3
! ! ! ! 				endif
! ! ! ! 
! ! ! ! 				Tkp = itk1 + sk1
! ! ! ! 				Tk1p= itk2 + sk2
! ! ! ! 				Tk2m= itk1 - sk1
! ! ! ! 				Tk3m= itk2 - sk2
! ! ! ! 			
! ! ! ! 				pkm1= 2.D0*( (itk1-1) - Tkp )*( (itk1-1) - Tkp )
! ! ! ! 				pk  = 2.D0*(  itk1    - Tk1p)*(  itk1    - Tk1p)
! ! ! ! 				pkp1= 2.D0*(  itk2    - Tk2m)*(  itk2    - Tk2m)
! ! ! ! 				pkp2= 2.D0*( (itk2+1) - Tk3m)*( (itk2+1) - Tk3m)
! ! ! ! 	
! ! ! ! 				t0= (t(it) + 2.d0 - Tkp )/( itk1-1 - Tkp )
! ! ! ! 				t1= (t(it) + 2.d0 - Tk1p)/( itk1   - Tk1p)
! ! ! ! 				t2= (t(it) + 2.d0 - Tk2m)/( itk2   - Tk2m)
! ! ! ! 				t3= (t(it) + 2.d0 - Tk3m)/( itk2+1 - Tk3m)
! ! ! ! 
! ! ! ! 				if(t(it)+2.d0 .gt. Tkp) then
! ! ! ! 					A0=0.d0
! ! ! ! 				else
! ! ! ! 					A0= fpt(pkm1, t0)
! ! ! ! 				endif
! ! ! ! 					A1=  fpt(pk  , t1)
! ! ! ! 					A2=  fpt(pkp1, t2)
! ! ! ! 		
! ! ! ! 				if(t(it)+2.d0 .lt. Tk3m) then
! ! ! ! 					A3=0.d0
! ! ! ! 				else
! ! ! ! 					A3= fpt(pkp2, t3)
! ! ! ! 				endif
! ! ! ! 
! ! ! ! 				xtaux=(A0*xk0 + A1*xk1 + A2*xk2+ A3*xk3)/(A0+A1+A2+A3)
! ! ! ! 				ytaux=(A0*yk0 + A1*yk1 + A2*yk2+ A3*yk3)/(A0+A1+A2+A3)
! ! ! ! 
! ! ! ! 				PCgamaSx(ip,PCnnp(ip)) = PCgamaSx(ip,PCnnp(ip))+
! ! ! !      &					(xtaux - xt(it))/ deltaS
! ! ! ! 				PCgamaSy(ip,PCnnp(ip)) = PCgamaSy(ip,PCnnp(ip))+
! ! ! !      &					(ytaux - yt(it))/ deltaS
! ! ! ! 			enddo
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
		endif
	
	end do

! 	if(idadossensib .eq.1)then
! 		write(*,*)'xsplxy - PCGama'
! 		do i=1,npc
! 			write(*,*)(PCGama(i,j),j=1,PCnnp(i))
! 		enddo
! 	endif
	
	return
	end



      SUBROUTINE DETCC(FI,DFI,VCC,FRES,NEPE,I_ELET,IC,IDCASO,NFIX,KODE)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      INCLUDE 'param_dim.inc'

      DIMENSION FI(NX),DFI(NX),KODE(NX)

      INTEGER FRES
      INTEGER IDCASO(NCasosX,3),I_ELET(16,NELEX)

C	Introduz todas as CC como fluxo nulo
C       no contorno externo (NFIX)
      DO I=1,NFIX
         KODE(I)=1
         FI(I)=0.0
      END DO

C     Atribui aos nos referentes aos eletrodos os valores das 
C     condicoes de contorno e IMPRIME

C     Eletrodos caso de carregamento em questao
c        IEL1= IDCASO(IC,1)
c        IEL2= IDCASO(IC,2)

      WRITE(FRES,1000)
 1000 FORMAT(//2X,'CONDICOES DE CONTORNO'//,
     &         5X,'ELETRODO',3X,'NO',6X,'CODIGO',5X,
     *            'VALORES PRESCRITOS')
      SVCC= VCC
      do iel =1,2
         SVCC= -1.D0*SVCC
         DO I=1,NEPE
            ino= I_ELET(idcaso(ic,iel),I)
            FI(ino)= SVCC
            WRITE(FRES,950) idcaso(ic,iel),ino, KODE(ino),FI(ino)
 950        FORMAT(I5,5X,5X,I3,8X,I1,8X,E14.7)
         ENDDO
      enddo

c     Identifica o eletrodo de referencia.
        IREF= IDCASO(IC,3)

C     Calcula o no que sera tomado como referencia (Potencial Zero)
C     (na metade do eletrodo de referencia)

      noref= I_ELET(IREF,1+NEPE/2)
      KODE(noref) = 0
      FI(noref)   = 0.d0

      WRITE(FRES,950) IREF,noref, KODE(noref),FI(noref)

      RETURN
      END




      SUBROUTINE FMAT(G,H,GF,HF,X,Y,XM,YM,FI,DFI,B,alfa,KODE,INC,NG,
     &                nfix,NN,N_eq,iaval)
C
C     Calculo das sub-matrizes G e H e formacao do 
C     sistema de equacoes Ax=b
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      INCLUDE 'param_dim.inc'
      INTEGER       FDAT,FDCC,FDEX,FRES,FROT
      COMMON /narq/ FDAT,FDCC,FDEX,FRES,FROT

      COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA


      DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max),NUMPCI(ninc_max)
      integer VARIA(ninc_max)
      DIMENSION PGLIM(nparg_max,3)

      DIMENSION GF(NFIX,NFIX),HF(NFIX,NFIX)
      DIMENSION G(N_eq,N_eq),H(N_eq,N_eq)
      DIMENSION X(NX),Y(NX),XM(NX),YM(NX),FI(NX),DFI(NX)
      DIMENSION INC(NX,2),KODE(NX)


C Numero total de elementos
      NE= NN
C Numero do n� movel inicial
      NINI= NFIX+1
         
C     Calculo da Coordenada dos Pontos de Colocacao
C     (Pontos medios dos elementos.)
C
      IF(IAVAL.EQ.1) THEN
         DO  I= 1,NE
            XM(I)=0.5*( X(INC(I,1)) + X(INC(I,2)) )
            YM(I)=0.5*( Y(INC(I,1)) + Y(INC(I,2)) )
         END DO
      ELSE
         DO  I= NINI,NE
            XM(I)=0.5*( X(INC(I,1)) + X(INC(I,2)) )
            YM(I)=0.5*( Y(INC(I,1)) + Y(INC(I,2)) )
         END DO
      ENDIF        

C     Zera as matrizes em funcao da variacao do numero de equacoes 
C     e do solver direto(fill-in)
      do i=1,N_eq
         dfi(i)= 0.d0
         do j=1,N_eq
            h(i,j)=  0.d0
            g(i,j)=  0.d0
         enddo
      enddo

C     
C CALCULO DAS MATRIZES G E H
C

C     Flag para indicar que nao se trata de calculo de valores em pontos internos.
      IPTI=0

      IF (IAVAL.EQ.1) THEN

C     Loop sobre os pontos de colocacao
         DO  I=1,NFIX
            DFI(I)=0.
C     Loop sobre os elementos
            DO J=1,NFIX
               NI= INC(J,1)
               NF= INC(J,2)
               IF(I.EQ.J) THEN
                  CALL INLO(X(NI),Y(NI),X(NF),Y(NF),G(I,J),H(I,J),
     &                      BJ,alfa)
               ELSE
                  CALL INTE(XM(I),YM(I),X(NI),Y(NI),X(NF),Y(NF),H(I,J),
     &                 dh1,dh2,G(I,J),dg1,dg2,B,BJ,db1,db2,alfa,IPTI,NG)
               END IF
               DFI(I)= DFI(I) - BJ
            ENDDO
         ENDDO

C     Armazena a parte fixa das matrizes
         do i=1,nfix
            do j=1,nfix
               gf(i,j)= g(i,j)
               hf(i,j)= h(i,j)
            enddo
         enddo

      ELSE

C     COPIA AS PARTES FIXAS DAS MATRIZES 
C   *******  Dispens�vel no caso de solu��o iterativa *****
         do i=1,nfix
            do j=1,nfix
               g(i,j)= gf(i,j)
               h(i,j)= hf(i,j)
            enddo
         enddo
         
      ENDIF


C     Loop sobre os pontos de colocacao FIXOS
      DO  I=1,NFIX
C     Loop sobre os elementos M�VEIS
         DO J=NINI,NN
            NI= INC(J,1)
            NF= INC(J,2)
            IF(I.EQ.J) THEN
               CALL INLO(X(NI),Y(NI),X(NF),Y(NF),G(I,J),H(I,J),BJ,alfa)
            ELSE
               CALL INTE(XM(I),YM(I),X(NI),Y(NI),X(NF),Y(NF),H(I,J),
     &              dh1,dh2,G(I,J),dg1,dg2,B,BJ,db1,db2,alfa,IPTI,NG)
            END IF
            DFI(I)= DFI(I) - BJ
         ENDDO
      ENDDO




C     Loop sobre os pontos de colocacao MOVEIS
      
      DO  I=NINI,NN
C     Loop sobre os elementos FIXOS E MOVEIS
         DO J=1,NN
            NI= INC(J,1)
            NF= INC(J,2)
            IF(I.EQ.J) THEN
               CALL INLO(X(NI),Y(NI),X(NF),Y(NF),G(I,J),H(I,J),BJ,alfa)
            ELSE
               CALL INTE(XM(I),YM(I),X(NI),Y(NI),X(NF),Y(NF),H(I,J),
     &              dh1,dh2,G(I,J),dg1,dg2,B,BJ,db1,db2,alfa,IPTI,NG)
            END IF
            DFI(I)= DFI(I) - BJ
         ENDDO
      ENDDO

      
      
      
c     Introdu��o das condi��es de compatibilidade nas inclusoes

C     Inicializa numero de nos acumulado.
      nna = NN
      nna0= NFIX
      
      do iinc= 1,nincl
         
C     Numero de nos (elementos) da inclusao
         nni= nninc(iinc)
         
C Introduz a parte de G referente a interface

         do j=1,nni
            
C           Linhas at� NN
            do i = 1,NN
               h(i,nna+j)= -g(i,nna0+j)
            end do

C           Linhas entre NN+1 e N_Eq

            do i = 1,nni
c              linhas de G afetadas da rela��o entre condutividades

               h(nna+i,nna+j)= c0*g(nna0+i,nna0+j)/cinc(iinc)

c              Poderia ser obtido, alternativamente, a partir dos 
c              coeficientes de G j� introduzidos em H como:

c              h(nna+i,nna+j)= -c0*h(nna0+i,nna+j)/cinc(iinc)


C              Para o elemento constante:linhas de H com sinal invertido, 
C              a menos da diagonal

               if (i.eq.j) then
                  h(nna+i,nna0+j)=  h(nna0+i,nna0+j)
               else
                  h(nna+i,nna0+j)= -h(nna0+i,nna0+j)               
               endif
               
            end do
         end do


C Numero de nos acumulado ate a inclusao anterior

         nna = nna  + nni
         nna0= nna0 + nni

      enddo

C
C Inicio da INTRODUCAO DAS CC: troca das colunas (e sinais) 
C                              se u for prescrito (KODE=0).
C
      DO J=1,NFIX
         IF(KODE(J).EQ.0)THEN
            DO I=1,NN
               CH=G(I,J)
               G(I,J)=-H(I,J)
               H(I,J)=-CH
            ENDDO
         END IF
      ENDDO

C
C     DFI originalmente contem os valores do termo forcado
C     Depois  da solucao contera os valores das incognitas no contorno.
C
      DO  I=1,NN
         DO  J=1,NFIX
            DFI(I)=DFI(I)+G(I,J)*FI(J)
         ENDDO
      ENDDO


!       r=0.
!       do i=1,N_eq
!         do j=1,N_eq
!             r = r + H(i,j)**2
!         enddo
!       enddo
! 
!         write(*,*)'Apos montagem |H| = ',r
! 
! 	write(*,*)'N_Eq',N_eq
! 	open(2,file='axb.m')
! 	write(2,*)' A = ['
! 	do i=1,N_eq
! 		write(2,11)(H(i,j),j=1,N_eq),'; '
! 	enddo
! 	write(2,*)' ]'
! 
! 	write(2,*)' y = ['
! 	do i=1,N_eq
! 		write(2,11)DFI(i)
! 	enddo
! 	write(2,*)' ]'
! 
! 	close(2)
!  11     format(520f24.12,a1)

      RETURN
      END




      SUBROUTINE INTE(XP,YP,XI,YI,XF,YF,HJ,DHJ1,DHJ2,
     &                GJ,DGJ1,DGJ2,B,BJ,DBJ1,DBJ2,alfa,IPTI,NG)
C
C ESTA SUBROTINA CALCULA O VALOR DOS ELEMENTOS FORA DA DIAGONAL DA
C MATRIZ H E G ATRAVES DE INTEGRACAO NUMERICA AO LONGO DOS ELEMENTOS
C DO CONTORNO.
C

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      COMMON /PG/ GI(8,2),OME(8,2)

      IF(NG.EQ.4)THEN
         ING=1
      ELSE IF(NG.EQ.8) THEN
         ING=2
      ELSE
         STOP'Numero incorreto de Pontos de Gauss'
      END IF
      
      AX= 0.5*(XF-XI)
      AY= 0.5*(YF-YI)

C     Calculo de L/2
      XL2=DSQRT(AX*AX + AY*AY)

C     Calcula (r*n)L/2 : (r escalar n)* L/2 
C     (constante para o elemento constante)
      RNXL2= (XI-XP)*AY - (YI-YP)*AX 


      Xmed= 0.5*(XF+XI)
      Ymed= 0.5*(YF+YI)

      GJ=   0.D0
      DGJ1= 0.D0
      DGJ2= 0.D0

      HJ=   0.D0
      DHJ1= 0.D0
      DHJ2= 0.D0

      BJ=   0.D0
      DBJ1= 0.D0
      DBJ2= 0.D0

      DO 40 I=1,NG

C     Coordenadas do ponto de integracao: 
         XCO= Xmed + AX*GI(I,ING)
         YCO= Ymed + AY*GI(I,ING)

C     Componentes do Raio
         Rx= XCO - XP
         Ry= YCO - YP
         Raio2= Rx*Rx + Ry*Ry
         Raio=  DSQRT(Raio2)
c         RNXL2= (Rx*AY - Ry*AX)

C     Integracao de G e H

         GJ= GJ + (alfa -  DLOG(Raio))      * OME(I,ING)
         HJ= HJ - (RNXL2 / Raio2 ) * OME(I,ING)

         BJ= BJ + (0.5 - DLOG(Raio))*OME(I,ING)

C     Calculo das derivadas em 
         IF (IPTI.EQ.1) THEN
            
            DGJ1= DGJ1 + Rx* OME(I,ING)/Raio2
            DGJ2= DGJ2 + Ry* OME(I,ING)/Raio2
            
            DHJ1= DHJ1 + ( AY - 2.d0*RNXL2*Rx/Raio2)* OME(I,ING)/Raio2
            DHJ2= DHJ2 + (-AX - 2.d0*RNXL2*Ry/Raio2)* OME(I,ING)/Raio2

	   DBJ1= DBJ1 + (RNXL2*Rx/Raio2-(0.5d0-DLOG(Raio))*AY)*OME(I,ING)
	   DBJ2= DBJ2 + (RNXL2*Ry/Raio2+(0.5d0-DLOG(Raio))*AX)*OME(I,ING)
	END IF


 40   CONTINUE

C        Multiplica lj/2 a Gij
      GJ= GJ * XL2

C     Calculo das derivadas em 
      IF (IPTI.EQ.1) THEN
         DGJ1= DGJ1 * XL2
         DGJ2= DGJ2 * XL2
      END IF

c        Multiplica B.r.n por bj
      BJ=   0.5D0 * B * RNXL2 * BJ
      DBJ1= 0.5D0 * B *         DBJ1
      DBJ2= 0.5D0 * B *         DBJ2

      RETURN
      END



      SUBROUTINE INLO(XI,YI,XF,YF,GJ,HJ,BJ,alfa)
C
C     Calculo dos valores dos elementos da diagonal de G e H.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (PI= 3.1415926535897932385D0)


      AXs2= 0.5d0*(XF-XI)
      AYs2= 0.5d0*(YF-YI)
      XLs2= DSQRT( AXs2*AXs2 + AYs2*AYs2 )
      GJ= 2.0d0 * XLs2 * (1.d0 - DLOG(XLs2) + alfa  )
      HJ= PI
      BJ= 0.0

      RETURN
      END




      SUBROUTINE SLNPDM(A,B,etol,d,IT,itri,N,nc2m)
C
C Solucao do sistema linear de equacoes pelo metodo de eliminacao de
C GAUSS, tendo a opcao de resolver sistemas ja triangularizados
C quando (ITRI.NE.0).
C
C
C A: Matriz do sistema.
C B: Matriz dos termos independentes.
C   Obs.: Originalmente B contem os coeficientes independentes e depois
C         de resolvido ela contem os valores das incognitas do sistema.
C IT: Vetor que acumula as linhas trocadas. E usado na resolucao de
C     sistemas ja triangularizados
C
C N: Numero atual de incognitas.
C


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DIMENSION A(N,N),B(N,NC2m),IT(N)


!         do j=1,nc2m
!             r=0.
!             do i=1,N
!                     r = r + b(i,j)**2
!             enddo
!         
!             write(*,*)j,'Antes solve slnpdm |b| = ',r
!         enddo

      N1=N-1

!       IF (ITRI.EQ.0) THEN
!          NC2M= 1
!       ELSE
!          NC2M= 2*NNP
!       END IF

      DO 10 K=1,N1
         K1=K+1
         C=A(K,K)
         IF (ITRI.EQ.0) THEN
            IT(K)= 0
         END IF
         IF(ABS(C).LE.etol) THEN
           DO 20 J= K1,N
C
C TENTA TROCAR LINHAS PARA OBTER COEFICIENTES DA DIAGONAL NAO NULOS.
C
              IF(ABS(A(J,K)).GT.etol) THEN
                 IT(K)= J
                 DO 30 L=K,N
                    C= A(K,L)
                    A(K,L)= A(J,L)
                    A(J,L)= C
 30              CONTINUE

                 DO 40 L=1,NC2M
                    C= B(K,L)
                    B(K,L)= B(J,L)
                    B(J,L)= C
 40              CONTINUE

                 C= A(K,K)
              ELSE

                 IF(K1.EQ.N) THEN
                   !  write(*,1001) K
 1001              FORMAT('**** SINGULARIDADE NA LINHA',I5)
                   stop
                   D= 0.
                   GO TO 999
                 END IF
              END IF
 20        CONTINUE
         END IF

C
C DIVIDE A LINHA PELO COEFICIENTE DA DIAGONAL
C
         IF (ITRI.EQ.0) THEN
            DO 50 J=K1,N
               A(K,J)=A(K,J)/C
 50         CONTINUE
         ELSE IF (IT(K).NE.0) THEN
            DO 60 J=1,NC2M
               CH= B(K,J)
               B(K,J)= B(IT(K),J)
               B(IT(K),J)= CH
 60         CONTINUE
         END IF

         DO 70 J=1,NC2M
            B(K,J)=B(K,J)/C
 70      CONTINUE

C
C ELIMINACAO DAS INCOGNITAS X(I) DA COLUNA K
C
         DO 80 I=K1,N
            C=A(I,K)
            IF (ITRI.EQ.0) THEN
               DO 90 J=K1,N
                  A(I,J)=A(I,J)-C*A(K,J)
 90            CONTINUE
            END IF

            DO 100 J=1,NC2M
               B(I,J)=B(I,J)-C*B(K,J)
 100        CONTINUE
 80      CONTINUE

 10   CONTINUE

C
C CALCULO DA ULTIMA INCOGNITA
C
      IF(ABS(A(N,N)).LT.etol) THEN
        !  write(*,1001)K
        D= 0.
        GOTO 999
      END IF

      DO 110 J=1,NC2M
         B(N,J)=B(N,J)/A(N,N)
 110  CONTINUE

C
C APLICACAO DO METODO DE RETROSUBSTITUICAO PARA CALCULO DAS INCOGNITAS
C RESTANTES
C
      DO 120 L=1,N1
         K=N-L
         K1=K+1
         DO 130 JJ=1,NC2M
            DO 140 J=K1,N
               B(K,JJ)=B(K,JJ)-A(K,J)*B(J,JJ)
 140        CONTINUE
 130     CONTINUE
 120  CONTINUE

C
C CALCULO DO VALOR DO DETERMINANTE
C
      D=1.
      DO 150 I=1,N
         D=D*A(I,I)
 150  CONTINUE

c      !  write(*,*) 'Determinante em SLNPDM',d

!         do j=1,nc2m
!             r=0.
!             do i=1,N
!                     r = r + b(i,j)**2
!             enddo
!         
!             write(*,*)j,'Apos solve slnpdm |b| = ',r
!         enddo


 999  RETURN
      END




      SUBROUTINE OUTPT(XM,YM,FI,DFI,KODE,NFIX,NN,N_eq)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'param_dim.inc'


      INTEGER       FDAT,FDCC,FDEX,FRES,FROT
      COMMON /narq/ FDAT,FDCC,FDEX,FRES,FROT
      COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA

      DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max),NUMPCI(ninc_max)
      integer VARIA(ninc_max)
      DIMENSION PGLIM(nparg_max,3)


      DIMENSION XM(NX),YM(NX),FI(NX),DFI(NX)
      DIMENSION KODE(NX)

      ZERO= 0.0
C
C REORDENAR OS VETORES FI E DFI PARA COLOCAR TODOS OS VALORES DO PO-
C TENCIAL EM FI E TODOS OS VALORES DA DERIVADA EM DFI
C
      DO 5 I=1,NFIX
         IF(KODE(I).GT.0) THEN 
           CH=FI(I)
           FI(I)=DFI(I)
           DFI(I)=CH
         END IF
 5    CONTINUE

C     Copia os valores calculados do 
C     potencial nas interfaces para FI

         nna = nn
         nna0= nfix

      do iinc =1,nincl
         nni = nninc(iinc)
         
         do i=1,nni
            fi(nna0+i)= dfi(nna0+i)
            dfi(nna0+i)= dfi(nna+i)
        enddo

         nna = nna  + nni
         nna0= nna0 + nni
      enddo

      WRITE(FRES,100)
  100 FORMAT(' ',80('*')//1X,'RESULTADOS'//,
     *       5X,'NOS DO CONTORNO'//3X,'No',6x,'X',
     *       13X,'Y',9X,'POTENCIAL (u)',7X,'D(u)/Dn'/)

      DO 10 I=1,NN
         WRITE(FRES,200)I,XM(I),YM(I),FI(I),DFI(I)
  200    FORMAT(I5,2(5X,E10.4),2(5X,E12.5))
   10 CONTINUE


      WRITE(FRES,600)
  600 FORMAT(' ',80('*'))

      icaso=1

C     IMPRESSAO DOS RESULTADOS PARA A OTIMIZACAO
      
c      WRITE(FROT,*) (FI(IP_OT(icaso,I)),I=1,14)

      RETURN
      END





 
      
!       subroutine centroide(x,y,N,NFIX,xc,yc,area)
! 
!       include 'param_dim.inc'
! 
!       integer N,NFIX
!       double precision x(NX),y(NX),area,area_parc,gxi,gyi,gx_sum
!       double precision xc,yc,xi,yi,xf,yf
!       
!       area=0
!       gx_sum=0
!       gy_sum=0
! 
!       
!       do i=NFIX+1,N-1
!          xi=x(i)
!          yi=y(i)
!          xf=x(i+1)
!          yf=y(i+1)
! 
!          area = area + area_parc(xi,yi,xf,yf)
!          gx_sum = gx_sum + gxi(xi,yi,xf,yf)
!          gy_sum = gy_sum + gyi(xi,yi,xf,yf)
!          
!       enddo
! 
!       xi=x(N)
!       yi=y(N)
!       xf=x(NFIX+1)
!       yf=y(NFIX+1)
! 
!       area =(area + area_parc(xi,yi,xf,yf))
!       gx_sum = gx_sum + gxi(xi,yi,xf,yf)
!       gy_sum = gy_sum + gyi(xi,yi,xf,yf)
!      
!       xc = gx_sum/area
!       yc = gy_sum/area
!      
!       
!       return
!       end



      subroutine plota_cont(p,fret,iter,numpar,numpc)

      INCLUDE 'param_dim.inc'
      save pz,xmin,xmax,ymin,ymax,iter_ant,ref_inic

      common/plota/ x,y,inc,N,NFIX,i_elet,nepe
      COMMON /cFDAT/ FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES

! 	common/plotcav/NNINC,NINCL
	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA
      common /NPAR/ NPAR
      
      CHARACTER*11 FILE_B,FILE_BD
      CHARACTER*15 FILE_DAT,FILE_DEX,FILE_OPT,FILE_RES
      character*18 FILE_PLT,FILE_PLD,FILE_JPG,FILE_PNG,INICIAL,EXATO
      CHARACTER*16 FILE_PLDE,FILE_PLDF 
      character*14 FILE_P
      character*1 c_num1,barra,ip,pl,sc
      character*2 c_num2
      character*26 ref
      character*26 ref_inic
      character*7 xp,yp,xch,ych
      character*9 resid,resid_inic
      character*60 label
          
      integer NNINC(ninc_max),NINCL
	double precision CINC(ninc_max),C0
	integer NUMPI(ninc_max), NUMPCI(ninc_max)
	integer VARIA(ninc_max)
      INTEGER FPLT,FPLD,FPLDE,FPLDF

      DOUBLE PRECISION X(NX),Y(NX),compr,xc,yc,area,Ju,Jv,Juv
      double precision xcinc(ninc_max),ycinc(ninc_max)
      double precision areainc(ninc_max),comprinc(ninc_max)
      double precision Juinc(ninc_max),Jvinc(ninc_max),Juvinc(ninc_max)	

      INTEGER INC(NX,2),I_ELET(16,NELEX)
      double precision p(nparg_max), pz(npar_max)

   
      FPLT= 20
      FPLD= 21
      FPLDE= 22
      FPLDF= 23

      data iplota/-1/
      iplota=iplota+1

      barra= achar(92)
      pl= achar(39)
!       write(resid,'(f7.4)')fret
      write(resid,'(e9.2)')fret

      FILE_B= FILE_OPT

      if (iter.EQ.0)then
         iter_ant=-1
!          write(resid_inic,'(f7.4)')fret
         write(resid_inic,'(e9.2)')fret
         ref_inic='Inicial - Resid.='//resid_inic 
      endif

        
      if (iter.EQ.iter_ant)then ! iteracao simples ou composta (powell)
         sc=' '
      else
         sc=' '
      endif

      iter_ant=iter

      

      if (iter.LE.9) then
            write (c_num1,'(I1)') iter
            ref='Iter.0'//c_num1//' - Resid.='//resid
         else
            write (c_num2,'(I2)') iter
            ref='Iter.'//c_num2//' - Resid.='//resid
      endif
 
       if (iplota.LE.9) then
          write (c_num1,'(I1)') iplota
          FILE_PLT= FILE_B//'-0'//c_num1//'.plt'
       else
          write (c_num2,'(I2)') iplota
          FILE_PLT= FILE_B//'-'//c_num2//'.plt'
       endif
 

  
      FILE_P= FILE_PLT
      FILE_JPG= FILE_P//'.jpg'
      FILE_PLD= FILE_P//'.pld'
      FILE_PNG= FILE_P//'.png'

      FILE_BD=FILE_DEX
      EXATO=FILE_BD//'-00.pld'

      INICIAL=FILE_B//'-00.pld'

      FILE_PLDE= FILE_B//'.plde'
      FILE_PLDF= FILE_B//'.pldf'


      OPEN (FPLT,FILE=FILE_PLT,STATUS='UNKNOWN')

      OPEN (FPLD,FILE=FILE_PLD,STATUS='UNKNOWN')
      

!        NFIX=N-NPAR

       if (iter.EQ.0) then
		do i=1,NFIX
	
			if(i.eq.1) then
				xmin= x(i)
				xmax= x(i)
				ymin= y(i)
				ymax= y(i)
			else if(i.le.N) then
				if(xmin.gt.x(i)) xmin= x(i)
				if(ymin.gt.y(i)) ymin= y(i)
				if(xmax.lt.x(i)) xmax= x(i)
				if(ymax.lt.y(i)) ymax= y(i)
			end if
		
		end do

		dx = xmax - xmin
		dy = ymax - ymin
	
		if (dx.gt.dy)then
			ymin = ymin-(dx-dy)/2.
			ymax = ymax+(dx-dy)/2.
		else if (dy.gt.dx)then
			xmin = xmin-(dy-dx)/2.
			xmax = xmax+(dy-dx)/2.
		endif

		dx = xmax - xmin
		xmin = xmin - .15*dx
		xmax = xmax + .15*dx
		ymin = ymin - .15*dx
		ymax = ymax + .15*dx

              
          OPEN (FPLDE,FILE=FILE_PLDE,STATUS='UNKNOWN')
          OPEN (FPLDF,FILE=FILE_PLDF,STATUS='UNKNOWN')
          write(FPLDF,*)'# PARTE FIXA DO CONTORNO'
          do i=1,NFIX
             write(FPLDF,*) X(inc(i,1)),Y(inc(i,1))
             write(FPLDF,*) X(inc(i,2)),Y(inc(i,2))
          end do
       
          do ielet=1,16
             do ielem=1,nepe
                i= i_elet(ielet,ielem)
                write(FPLDE,*) X(inc(i,1)),Y(inc(i,1))
                write(FPLDE,*) X(inc(i,2)),Y(inc(i,2))
             enddo
             write(FPLDE,*)
          end do
          close (FPLDE)
          close (FPLDF)
       
       endif

      
       
       write (FPLD,*)'# PARTE MOVEL DO CONTORNO'
       do i=NFIX+1,N
           write(FPLD,*) X(inc(i,1)),Y(inc(i,1))
           write(FPLD,*) X(inc(i,2)),Y(inc(i,2))
  	   write(FPLD,*)
       end do



       write(FPLD,*)
       write (FPLD,*)'# Pontos de controle'
       
! 	write(*,*)'nincl',nincl
	ipc = 0
	do j=1,nincl
		do i=1,numpci(j)
			write(FPLD,*) p(3*ipc + 3*i-2),p(3*ipc + 3*i-1),p(3*ipc + 3*i)
		enddo
		write(FPLD,*)p(3*ipc + 1),p(3*ipc + 2),p(3*ipc + 3) 
		ipc = ipc + numpci(j)
		write(FPLD,*)
	enddo
       

	nfimant=NFIX
	do i=1,nincl
		nfimatual=nfimant+NNINC(i)
      		call propgeom(x,y,nfimatual,nfimant,compr,
     &				area,xc,yc,Ju,Jv,Juv)
!        		call centroide(x,y,nfimatual,nfimant,xc,yc,area)
		xcinc(i)=xc
		ycinc(i)=yc
		comprinc(i)=compr
		areainc(i)=abs(area)
		Juinc(i) = Ju
		Jvinc(i) = Jv
		Juvinc(i) = Juv

		nfimant=nfimatual
	enddo

        write(FPLD,*)'#   Centroide (x,y)   '

	do i=1,nincl
 		write (FPLD,*) xcinc(i),ycinc(i)
 		write (FPLD,*) 
	enddo
       write(xch,'(f7.3)') xc-.2
       write(ych,'(f7.3)') yc

       write(FPLD,*)
       write(FPLD,*)'#  xc, yc, area, perimetro, Ju, Jv, Juv:'
	do i=1,nincl
 		write (FPLD,1070) ' #  ',i,xcinc(i),ycinc(i),
     &		areainc(i),comprinc(i),Juinc(i),Jvinc(i),Juvinc(i)
	enddo
       close (FPLD)
       
        open(1,file='pgeom.tmp') !so para uma inclusao
	write (1,*) xcinc(1),ycinc(1),sqrt(xcinc(1)**2+ycinc(1)**2),
     &		areainc(1),comprinc(1),Juinc(1),Jvinc(1)
        close(1)

 1070  format(a4,i5,7f20.4)


       if (iter.EQ.0) then
       
	  open(1,file='fobj0.tmp')
	  write(1,*)resid
	  close(1)       
          
          ref='Inicial - Resid.='//resid

c          label='set label '//pl//'G'//pl//' right front at '//
c     &                xch//','//ych//' tc lt 10'
c          write(FPLT,1040)label

          do i=1,numpc
             write(ip,'(I1)') i
             write(xp,'(f7.3)') p(3*i-2)+.2  
             write(yp,'(f7.3)') p(3*i-1)
             label= 'set label '//pl//ip//pl//
     &           ' left front at '//xp//','//yp//' tc lt 10'
             write(FPLT,1040)label
          enddo          
     
          write(FPLT,1050) FILE_B,
     &                  xmin,xmax,ymin,ymax,barra,
     &                  FILE_PLDF, barra, FILE_PLDE,barra,
     &                  EXATO,'Exato', barra,
     &                  FILE_PLD,ref,
     &                  './jpg/'//file_jpg,'./png/'//file_png


 1050     format('set key bottom', /,'set title',1h',a11,1h',/,
     &         'set size square',/,  
     &         'plot [',f10.2,':',f10.2,'][',f10.2,':',f10.2,']',a1,/,
     &            5x,1h',a16,1h',' notitle w linesp ,',a1,/,
     &            5x,1h',a16,1h',' notitle w linesp ,',a1,/,
     &  5x,1h',a18,1h',' t ',1h',a5,1h', ' w linesp lt 4 ps .4,',a1,/,
     &  5x,1h',a18,1h',' t ',1h',a26,1h',' w linesp lt 10 ps .4 ',/,
     &       'set terminal jpeg large size 600,600',/,
     &       'set output ',1h',a24,1h',/,'replot',/,
     &       '#set terminal png small size 400,400',/,
     &       '#set output ',1h',a24,1h',/,'replot')
       
       else
          
c          label='set label '//pl//'G'//pl//' right front at '//
c     &                   xch//','//ych//' tc lt 3'
c          write(FPLT,1040)label
	  open(1,file='fobjf.tmp')
	  write(1,*)resid
	  close(1)    
          do i=1,numpc
!              write(ip,'(I1)') i
!              write(xp,'(f7.3)') pz(3*i-2)+.2
!              write(yp,'(f7.3)') pz(3*i-1)
!              label= 'set label '//pl//ip//pl//
!      &              ' left front at '//xp//','//yp//' tc lt 10'
!              write(FPLT,1040)label
       
             write(ip,'(I1)') i
             write(xp,'(f7.3)') p(3*i-2)-.2
             write(yp,'(f7.3)') p(3*i-1)
             label= 'set label '//pl//ip//pl//
     &           ' right front at '//xp//','//yp//' tc lt 3'
             write(FPLT,1040)label
          enddo
 1040     format(a60)
    
          write(FPLT,1020) FILE_B,
     &                  xmin,xmax,ymin,ymax,barra,
     &                  FILE_PLDF, barra, FILE_PLDE,barra,
     &                  EXATO,'Exato', barra,
     &                  INICIAL,ref_inic, barra,
     &                  FILE_PLD,ref,
     &                  './jpg/'//file_jpg,'./png/'//file_png

 1020     format('set key bottom', /,'set title',1h',a11,1h',/,
     &         'set size square',/,  
     &         'plot [',f10.2,':',f10.2,'][',f10.2,':',f10.2,']',a1,/,
     &            5x,1h',a16,1h',' notitle w linesp ,',a1,/,
     &            5x,1h',a16,1h',' notitle w linesp ,',a1,/,
     &  5x,1h',a18,1h',' t ',1h',a5,1h', ' w linesp lt 4 ps .4,',a1,/,
     &  5x,1h',a18,1h',' t ',1h',a26,1h', ' w linesp lt 10 ps .4,',a1,/,
     &  5x,1h',a18,1h',' t ',1h',a26,1h',' w linesp lt 3 ps .4',/,
     &       'set terminal jpeg large size 600,600',/,
     &       'set output ',1h',a24,1h',/,'replot',/,
     &       '#set terminal png small size 400,400',/,
     &       '#set output ',1h',a24,1h',/,'replot')
       endif
     
       close(FPLT)

      call system('gnuplot '//FILE_PLT )


      return
      end








	SUBROUTINE FMATD(X,Y,XM,YM,DB,FI,DFI,alfa,KODE,INC,NG,NNP,NP,
     &			NFIX,NN,N_eq)

	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	include 'param_dim.inc'

	DIMENSION X(NX),Y(NX),XM(NX),YM(NX),FI(NX),DFI(NX),
     *          DH(NX,NX),DG(NX,NX),DB(N_eq,NNP*2),NP(NX)
	DIMENSION INC(NX,2),KODE(NX)
	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA

	DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max),NUMPCI(ninc_max)
	integer VARIA(ninc_max)
	DIMENSION PGLIM(nparg_max,3)

	DO 10 IP=1,NNP !loop dos nohs parametros


		NEPP= NP(IP) !numero do elemento posterior ao noh parametro
		NEAP= NEPP -1 !numero do elemento anterior ao noh parametro


		if(INC(NEAP,2).ne.NP(IP))then  ! procurar NEAP se o noh parametro for o 1o de cada spline
			do j=NFIX+1,NN
				if(INC(j,2).eq.NP(IP))then
					NEAP = j
				endif
			enddo
		endif



		do j=1,N_eq  !inicialização das matrizes DH e DG
			do i=1,N_eq
				DG(i,j)=0.
				DH(i,j)=0.
			enddo
		enddo



		DO 20 ID=1,2   ! derivada nas duas direcoes

			IC= (IP-1)*2 + ID !coluna da matriz DB referente ao noh parmetro IP, direcao ID

			DO 30 I=1,NN  ! loop sobre os pontos de colocacao (fixos e moveis)
				DO 40 J=1,NN ! loop sobre os elementos (fixos e moveis)
					INFA= 0 ! ponto de colocacao em elemento alterado INFA = 1
					INCA= 0 ! integracao em elemento alterado INCA = 1
					IF (I.EQ.NEPP) THEN  ! ponto de colocacao em elemento alterado pelo noh parametro

						INFA= 1
						IF (I.EQ.J) THEN
							CALL INLOD(X(inc(J,1)),Y(inc(J,1)),
     &								   X(inc(J,2)),Y(inc(J,2)),
     *								   DG(I,J),NEAP,ID,J,alfa)
							DH(I,J)= 0.
						ELSE IF ( J.EQ.NEAP) THEN
							INCA= 1
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *							DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)
						ELSE
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *  						DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)
                     				END IF

					ELSE IF (I.EQ.NEAP) THEN !ponto de colocacao em elemento alterado pelo noh parametro

                    				INFA= 1
						IF (I.EQ.J) THEN
							CALL INLOD(X(inc(J,1)),Y(inc(J,1)),
     &								   X(inc(J,2)),Y(inc(J,2)),
     *								DG(I,J),NEAP,ID,J,alfa)
							DH(I,J)= 0.
						ELSE IF ( J.EQ.NEPP) THEN
							INCA= 1
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *							DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)

						ELSE
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *							DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)
						END IF

					ELSE IF ((J.EQ.NEAP).OR.(J.EQ.NEPP)) THEN !ponto de colocacao em elemento NAO alterado e integracao sobre elementos alterados
						INCA= 1
						CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &								       X(inc(J,2)),Y(inc(J,2)),
     *  					DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)

					END IF

 40            			CONTINUE
 30         		CONTINUE

! compatibilidades nas inclusoes (analogo a montagem das matrizes H e G)

			nna = NN ! inicializa o numero de nos acumulado.
			nna0= NFIX

			do iinc= 1,nincl
				nni= nninc(iinc)  !Numero de nos (elementos) da inclusao

				do j=1,nni
					do i = 1,NN ! linha de DH ate NN
						DH(i,nna+j)= -DG(i,nna0+j)!copia as colunas de DG referentes ao fluxo nos elementos das inclusoes
					end do

					do i = 1,nni ! linhas de DH de NN ate N_eq
						DH(nna+i,nna+j)= c0*DG(nna0+i,nna0+j)/cinc(iinc) ! coloca em DH as linhas de DG afetadas da relacao entre condutividades
						if (i.eq.j) then
							DH(nna+i,nna0+j)=  DH(nna0+i,nna0+j) !DH eh sempre zero na diagonal
						else
							DH(nna+i,nna0+j)= -DH(nna0+i,nna0+j) !copia parte de DH com sinal invertido
						endif
					end do
				end do

				nna = nna  + nni ! atualiza o numero de nohs acumulados
				nna0= nna0 + nni
			enddo !fim do loop da relacao de compatibilidades


			DO 60 J=1,NFIX !Rearranjo das matrizes DH(elemntos referentes as incog ja calc DFI) e DG(valores prescritos FI)
				IF(KODE(J).EQ.0) THEN
					DO 70 I=1,NN !ate NN pq de NN ate N_eq 
						CH= DH(I,J)
						DH(I,J)= -DG(I,J)
						DG(I,J)= -CH
 70               			CONTINUE
				END IF
 60         		CONTINUE


			do i=1,N_eq
				DB(I,IC)= 0.
			enddo


			DO J=1,NFIX ! FI possui apenas NFIX linhas
				DO I=1,N_eq ! Montagem da coluna IC
					DB(I,IC)= DB(I,IC) + DG(I,J)*FI(J)
				enddo
			ENDDO
			DO J=1,N_eq ! DFI possui N_eq linhas
				DO I=1,N_eq ! Montagem da coluna IC
					DB(I,IC)= DB(I,IC) - DH(I,J)*DFI(J)
				enddo				
			ENDDO


 20     	CONTINUE !fim do loop das direcoes
 10   	CONTINUE !fim do loop dos nohs parametros 



      RETURN
      END

      SUBROUTINE INTED(XP,YP,X1,Y1,X2,Y2,DH,DG,NG,INCA,INFA,
     &			NEAP,ID,J,alfa)
C
C Esta subrotina calcula o valor dos elementos fora da diagonal da
C matriz dH e dG atraves de integracao numerica ao longo dos elementos
C do contorno.
C
C R= Distancia do ponto em consideracao ate os pontos de integracao nos
C elementos de contorno.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DOUBLE PRECISION L, NX, NY, JE,alfa

      DIMENSION XCO(8),YCO(8),GI(8),OME(8),DR(2),L(2)

      IF (NG.EQ.8) THEN

         GI(1)=0.96028986
         GI(2)=-GI(1)
         GI(3)=0.79666648
         GI(4)=-GI(3)
         GI(5)=0.52553210
         GI(6)=-GI(5)
         GI(7)=0.18343464
         GI(8)=-GI(7)

         OME(1)=0.10122854
         OME(2)=OME(1)
         OME(3)=0.22238103
         OME(4)=OME(3)
         OME(5)=0.31370665
         OME(6)=OME(5)
         OME(7)=0.36268378
         OME(8)=OME(7)

      ELSE

         GI(1)=0.86113631
         GI(2)=-GI(1)
         GI(3)=0.33998104
         GI(4)=-GI(3)

         OME(1)=0.34785485
         OME(2)=OME(1)
         OME(3)=0.65214515
         OME(4)=OME(3)

      END IF

      BX= (X2 + X1) / 2.
      BY= (Y2 + Y1) / 2.

      XLX= (X2 - X1)
      XLY= (Y2 - Y1)

      AX= XLX / 2.
      AY= XLY / 2.

      XL= DSQRT( XLX**2 + XLY**2 )

      JE= XL / 2.

      NX= XLY / XL
      NY= -XLX / XL

      DG=0.
      DH=0.

      DO 10 I=1,NG

         XCO(I)= AX*GI(I) + BX
         YCO(I)= AY*GI(I) + BY

         RX= XCO(I) - XP
         RY= YCO(I) - YP

         R2=  RX**2 + RY**2
         R= DSQRT( R2 )

         UST=  alfa + DLOG(1/R)
         QST=  -(RX*NX + RY*NY) / R2

C
C  Calculo de dR, dN e dJe
C
      DR(1)= 0.
      DR(2)= 0.

      DNX= 0.
      DNY= 0.
      DJE= 0.

      L(1)= XLX
      L(2)= XLY

      IF (INCA.EQ.1) THEN

         SIGN= 1.
         IF(J.EQ.NEAP) SIGN= -1.

         DR(ID)= (1 - SIGN*GI(I)) / 2.
         DJE= -SIGN*L(ID) / (2*XL)

         IF (ID.EQ.1) THEN
            DNX= SIGN*(XLX*NX) / (XL*XL)
            DNY= SIGN*( (XLX*NY)/XL + 1)/XL
         ELSE
            DNX= SIGN*( (XLY*NX)/XL - 1)/XL
            DNY= SIGN*(XLY*NY) / (XL*XL)
         END IF
      END IF

      IF (INFA.EQ.1) THEN
         DR(ID)= DR(ID) - 0.5
      END IF

C
C  Calculo de dU* e dQ*
C

      A= (DR(1)*NX + DR(2)*NY + RX*DNX + RY*DNY) / R2

      B= ( (RX*NX +RY*NY) * (RX*DR(1) + RY*DR(2)) ) / (R2*R2)

      DQST= 2*B - A
      DUST= -(RX*DR(1) + RY*DR(2)) / R2

C
C  Calculo de dG e dH
C

         DG= DG + (JE*DUST + UST*DJE) * OME(I)
         DH= DH + (JE*DQST + QST*DJE) * OME(I)

 10   CONTINUE

      RETURN
      END


      SUBROUTINE INLOD(X1,Y1,X2,Y2,DG,NEAP,ID,J,alfa)

C
C ESTA SUBROTINA CALCULA O VALOR DOS ELEMENTOS DA DIAGONAL DA MATRIZ dG
C

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DOUBLE PRECISION L

      DIMENSION L(2)

      L(1)= (X2 - X1)
      L(2)= (Y2 - Y1)

      XL= DSQRT( L(1)**2 + L(2)**2)

      SIGN= 1.
      IF (J.EQ.NEAP) SIGN= -1.

      DG= -SIGN*(alfa+ DLOG( 2/XL ))*L(ID)/XL

      RETURN
      END


      subroutine propgeom(x,y,N,NFIX,compr,area,xc,yc,Ju,Jv,Juv)

      include 'param_dim.inc'

      integer N,NFIX
      double precision x(N),y(N),area,area_parc,gxi,gyi
      double precision gx_sum,gy_sum,compr
      double precision Jx,Jy,Jx_parc,Jy_parc,Jxy,Jxy_parc,Ju,Jv,Juv
      double precision xc,yc,xc1,yc1,xi,yi,xf,yf,alfa

    
      area=0
      compr=0
      gx_sum=0
      gy_sum=0
      Jx = 0
      Jy = 0
      Jxy = 0
      
      do i=NFIX+1,N-1
         xi=x(i)
         yi=y(i)
         xf=x(i+1)
         yf=y(i+1)

         area = area + area_parc(xi,yi,xf,yf)
         compr= compr + sqrt((xf-xi)*(xf-xi)+(yf-yi)*(yf-yi))
         gx_sum = gx_sum + gxi(xi,yi,xf,yf)
         gy_sum = gy_sum + gyi(xi,yi,xf,yf)

      enddo

      xi=x(N)
      yi=y(N)
      xf=x(NFIX+1)
      yf=y(NFIX+1)

      area =-(area + area_parc(xi,yi,xf,yf))
      compr= compr + sqrt((xf-xi)*(xf-xi)+(yf-yi)*(yf-yi))

      gx_sum = -(gx_sum + gxi(xi,yi,xf,yf))
      gy_sum = -(gy_sum + gyi(xi,yi,xf,yf))

      xc = gx_sum/area
      yc = gy_sum/area
      

      
      do i=NFIX+1,N-1
         xi=x(i)  -xc1
         yi=y(i)  -yc1
         xf=x(i+1)-xc1
         yf=y(i+1)-yc1
         
         Jx = Jx + Jx_parc(xi,yi,xf,yf)
         Jy = Jy + Jy_parc(xi,yi,xf,yf)
         Jxy= Jxy+ Jxy_parc(xi,yi,xf,yf)
      enddo 
      
      xi=x(N)    -xc1
      yi=y(N)    -yc1
      xf=x(NFIX+1)-xc1
      yf=y(NFIX+1)-yc1
      
      Jx = -(Jx + Jx_parc(xi,yi,xf,yf))
      Jy = -(Jy + Jy_parc(xi,yi,xf,yf))  
      Jxy=-(Jxy+ Jxy_parc(xi,yi,xf,yf))

c     Calculo dos momentos de inercia segundo as direcoes...
c     principais 
c     lembrar que alfa=2*teta
c     Juv so para confirmar que da zero


        alfa=atan(-2*Jxy/(Jx-Jy))
	if(Jx .eq. Jy )then 
		alfa=0.
	endif


      Ju = (Jx+Jy)/2 + (Jx-Jy)*cos(alfa)/2 - Jxy*sin(alfa)
      Jv = (Jx+Jy)/2 - (Jx-Jy)*cos(alfa)/2 + Jxy*sin(alfa)
      Juv= (Jx-Jy)*sin(alfa)/2 + Jxy*cos(alfa)


      return
      end

      function area_parc(xi,yi,xf,yf)
      double precision area_parc,xi,yi,xf,yf
      
      area_parc=(yf* xi - xf*yi)/2
      
      return
      end

      
	function gxi(xi,yi,xf,yf)
	double precision gxi,xi,yi,xf,yf
	
		gxi=(xf*yf*xi - xf*xi*yi)/6 - (xf*xf*yf)/3
     &             +(yf*xi*xi)/6 - (xf*xf*yi)/6 + (xi*xi*yi)/3
	
	return
	end

	function gyi(xi,yi,xf,yf)
	double precision gyi,xi,yi,xf,yf
	
		gyi=(yf*xi*yi - xf*yf*yi)/6 + (xf*yf*yf)/3 - (xf*yi*yi)/6 
     &            + (yf*yf*xi)/6 - (xi*yi*yi)/3
	
	return
	end

	function Jx_parc(xi,yi,xf,yf)
	double precision Jx_parc,xi,yi,xf,yf,xc,yc
	
	Jx_parc= (xf*yf*yf*yf)/4 - (xf*yi*yi*yi)/12 + (yf*yf*yf*xi)/12
     &         - (xi*yi*yi*yi)/4 - (xf*yf*yi*yi)/12 - (xf*yf*yf*yi)/12
     &         + (yf*xi*yi*yi)/12 + (yf*yf*xi*yi)/12
	
	
	return
	end

      function Jy_parc(xi,yi,xf,yf)
      double precision Jy_parc,xi,yi,xf,yf,xc,yc
    
      Jy_parc= (yf*xi*xi*xi)/12 - (xf*xf*xf*yf)/4 - (xf*xf*xf*yi)/12 
     &       + (xi*xi*xi*yi)/4 + (xf*yf*xi*xi)/12 + (xf*xf*yf*xi)/12
     &       - (xf*xi*xi*yi)/12 - (xf*xf*xi*yi)/12
      
      return
      end

      function Jxy_parc(xi,yi,xf,yf)
      double precision Jxy_parc,xi,yi,xf,yf

      Jxy_parc= (xf*yf*yf*xi)/12 - (xf*xf*yf*yi)/12 - (xf*xi*yi*yi)/12 
     &        + (yf*xi*xi*yi)/12 + (xf*xf*yf*yf)/8  - (xf*xf*yi*yi)/24 
     &        + (yf*yf*xi*xi)/24 - (xi*xi*yi*yi)/8
      
      return
      end


	subroutine matvec(A,x,y,N)

	real*8 A(N,N),x(N),y(N)

	do i=1,N
		y(i) = 0.
	enddo

	do j=1,N
		do i=1,N
			y(i) = y(i) + A(i,j)*x(j)
		enddo
	enddo
	
	return
	end
   
